This notebook analyses the Part A data in terms of variable importance using a random forest model based on conditional inference trees and a conditional permutation variable importance algorithm.
# load packages
library(tidyr)
library(ggplot2)
library(party)
library(conflicted)
library(tidyverse)
library(openxlsx)
library(caret)
library(viridis)
library(cowplot)
# set package parameters
theme_set(theme_bw())
# plot colour scheme
mycolourlist = list(c(0, 102, 255), c(0, 204, 153), c(255, 0, 102), c(74, 111, 152), c(251, 164, 49), c(204, 153, 255), c(90, 192, 255), c(80, 245, 233), c(255, 90, 192), c(164, 201, 242), c(255, 254, 139), c(255, 243, 255))
mycolours = matrix()
for (ii in 1:length(mycolourlist)){
mycolours[ii] = rgb(mycolourlist[[ii]][1]/255,
mycolourlist[[ii]][2]/255,
mycolourlist[[ii]][3]/255)
}
# toggle to save plots
saveplots = TRUE
if (saveplots){
# set output plot directory
choose.files(caption="Just cancel this", filters=matrix(data=c(" ", " "), ncol=2)) # workaround for bug in RTerm choose.dir
outFigPath <- utils::choose.dir(caption="Select output folder to save plots '03 Experiment\\Experiment 1\\Analysis\\Plots'")
if (!dir.exists(file.path(outFigPath, "svg"))){dir.create(file.path(outFigPath, "svg"))}
if (!dir.exists(file.path(outFigPath, "pdf"))){dir.create(file.path(outFigPath, "pdf"))}
}
# toggle to save data
savedata = TRUE
if (savedata){
# set output plot directory
if (saveplots==FALSE){
choose.files(caption="Just cancel this", filters=matrix(data=c(" ", " "), ncol=2)) # workaround for bug in RTerm choose.dir
}
outDataPath <- utils::choose.dir(caption="Select output folder to save data '03 Experiment\\Experiment 1\\Analysis\\R'")
}
stimDataApath <- utils::choose.files(caption=r"(Select refmap_listest1_testdataA_ByStim.csv from 03 Experiment\Experiment 1\Analysis\PostProcess)",
filters=matrix(data=c("refmap_listest1_testdataA_ByStim.csv", "refmap_listest1_testdataA_ByStim.csv"), ncol=2))
stimNoticeDataApath <- utils::choose.files(caption=r"(Select refmap_listest1_testdataANoticeFilt_ByStim.csv from 03 Experiment\Experiment 1\Analysis\PostProcess)", filters=matrix(data=c("refmap_listest1_testdataANoticeFilt_ByStim.csv", "refmap_listest1_testdataANoticeFilt_ByStim.csv"), ncol=2))
stimDataA <- utils::read.csv(stimDataApath, header=TRUE)
stimNoticeDataA <- utils::read.csv(stimNoticeDataApath, header=TRUE)
colnames(stimDataA)[1] <- "Stimulus"
colnames(stimNoticeDataA)[1] <- "Stimulus"
# make response proportions into percentages
stimDataA[['HighAnnoyPc']] <- stimDataA[['HighAnnoyProp']]*100
stimDataA[['dHighAnnoyPc']] <- stimDataA[['dHighAnnoyProp']]*100
# make response proportions into percentages
stimNoticeDataA[['HighAnnoyPcFilt']] <- stimNoticeDataA[['HighAnnoyPropFilt']]*100
stimNoticeDataA[['dHighAnnoyPcFilt']] <- stimNoticeDataA[['dHighAnnoyPropFilt']]*100
stimNoticeDataA[['NoticedPcFilt']] <- stimNoticeDataA[['NoticedPropFilt']]*100
# function to encode categorical to ordinal numeric variables
encode_ordinal <- function(x, order=unique(x)) {
x <- as.numeric(factor(x, levels=order, exclude=NULL, order=TRUE))
x
}
# definition of ordinal variable levels
SNRCatsA <- c("No UAS", "-16", "-10", "-4", "2", "8")
UASLAeqCatsA <- c("No UAS", "42", "48", "54", "60")
The aggregated data by stimulus are assigned to a dataframe, relevant categorical variables are converted to ordinal, and then the variable subset of interest is selected, NA rows dropped (ie, the ‘no UAS’ stimuli, as the conditional variable importance algorithm cannot currently handle NA values, which are present in all the UAS dB metrics), and a formula assigned.
stimDataANum <- data.frame()
stimDataANum <- cbind(stimDataA[, 'Stimulus'],
stimDataA[, which(colnames(stimDataA)=="UASLAeq"):
which(colnames(stimDataA)=="SNRlevel")],
stimDataA[, "IntermitRatioMaxLR", drop=F],
stimDataA[, which(colnames(stimDataA)=="UASLAeqMaxLR"):
which(colnames(stimDataA)=="UASImpulsSHM05ExMaxLR")],
stimDataA[, which(colnames(stimDataA)=="LAeqLAF90diff"):
which(colnames(stimDataA)=="dImpulsSHM05ExMaxLR")],
stimDataA[, which(colnames(stimDataA)=="ValenceMedian"):
which(colnames(stimDataA)=="dHighAnnoyProp")],
stimDataA[, which(colnames(stimDataA)=="HighAnnoyPc"):
which(colnames(stimDataA)=="dHighAnnoyPc")])
# remove duplicated variables
stimDataANum <- subset(stimDataANum, select = -c(UASLAeqMaxLR, UASLAEMaxLR))
colnames(stimDataANum)[1] <- "Stimulus"
# make the categorical outcome variables numerical dummy factors
# stimDataANum[['SNRlevel']] <- encode_ordinal(stimDataA[['SNRlevel']], order=SNRCatsA)
# stimDataANum[['UASLAeq']] <- encode_ordinal(stimDataA[['UASLAeq']], order=UASLAeqCatsA)
# make the discrete ordinal outcome variables factors
stimDataANum[['ValenceMedian']] <- factor(stimDataANum[['ValenceMedian']], levels=c(1, 2, 3, 4, 5), order=TRUE)
stimDataANum[['ArousalMedian']] <- factor(stimDataANum[['ArousalMedian']], levels=c(1, 2, 3, 4, 5), order=TRUE)
stimDataANum[['AnnoyMedian']] <- factor(stimDataANum[['AnnoyMedian']], levels=c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), order=TRUE)
stimDataANum[['dValenceMedian']] <- factor(stimDataANum[['dValenceMedian']], levels=c(-4, -3, -2, -1, 0, 1, 2, 3, 4), order=TRUE)
stimDataANum[['dArousalMedian']] <- factor(stimDataANum[['dArousalMedian']], levels=c(-4, -3, -2, -1, 0, 1, 2, 3, 4), order=TRUE)
stimDataANum[['dAnnoyMedian']] <- factor(stimDataANum[['dAnnoyMedian']], levels=c(-10, -9, -8, -7, -6, -5, -4, -3, -2, -1,
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), order=TRUE)
stimNoticeDataANum <- data.frame()
stimNoticeDataANum <- cbind(stimNoticeDataA[, 'Stimulus'],
stimNoticeDataA[, which(colnames(stimNoticeDataA)=="UASLAeq"):
which(colnames(stimNoticeDataA)=="SNRlevel")],
stimNoticeDataA[, "IntermitRatioMaxLR", drop=F],
stimNoticeDataA[, which(colnames(stimNoticeDataA)=="UASLAeqMaxLR"):
which(colnames(stimNoticeDataA)=="UASImpulsSHM05ExMaxLR")],
stimNoticeDataA[, which(colnames(stimNoticeDataA)=="LAeqLAF90diff"):
which(colnames(stimNoticeDataA)=="dImpulsSHM05ExMaxLR")],
stimNoticeDataA[, which(colnames(stimNoticeDataA)=="NoticedTotalFilt"):
which(colnames(stimNoticeDataA)=="NoticedPcFilt")])
stimNoticeDataANum <- subset(stimNoticeDataANum, select = -c(UASLAeqMaxLR, UASLAEMaxLR))
colnames(stimNoticeDataANum)[1] <- "Stimulus"
# make the categorical outcome variables numerical dummy factors
# stimNoticeDataANum[['SNRlevel']] <- encode_ordinal(stimNoticeDataA[['SNRlevel']], order=SNRCatsA)
# stimNoticeDataANum[['UASLAeq']] <- encode_ordinal(stimNoticeDataA[['UASLAeq']], order=UASLAeqCatsA)
Write a function to train a conditional-inference random forest (crf) model on input data according to input formula, iterate over input random seeds, average error and variable importance metrics, and output metrics with plotted
multi_crfReg <- function(dataIn, iVars, dVar, seeds, ntree, mtry, permImpCondThres=0.95, minsplit=20, minbucket=7, nperm=1){
# initialise variables
crfOOBErrAll <- 0
crfOOBRMSE <- 0
crfOOBMAE <- 0
crfOOBErrR2 <- 0
crfMarPermImpVals <- 0
crfConPermImpVals <- 0
crfMarPermImpValsPerTree <- data.frame()
crfConPermImpValsPerTree <- data.frame()
for (iters in 1:length(seeds)){
# formula for regression
formVars <- reformulate(iVars, dVar)
# set random seed
set.seed(seeds[iters])
# train crf model
crfModel <- party::cforest(formVars, data=dataIn,
controls=party::cforest_unbiased(ntree=ntree,
mtry=mtry,
minsplit=minsplit,
minbucket=minbucket))
# get OOB predictions
crfModelOOB <- predict(crfModel, OOB=TRUE, type='response')
# get OOB error
crfModelOOBErr <- as.numeric(as.matrix(as.numeric(as.matrix(crfModelOOB))
- as.numeric(as.matrix(crfModel@data@env$response[[names(crfModel@data@env$response)]]))))
# OOB RMSE, error quartiles and Rsquared
crfOOBRMSE <- crfOOBRMSE + sqrt(mean(crfModelOOBErr^2))
crfOOBMAE <- crfOOBMAE + mean(abs(crfModelOOBErr))
crfOOBErrR2 <- crfOOBErrR2 + cor(as.numeric(as.matrix(crfModelOOB)),
as.numeric(as.matrix(crfModel@data@env$response[[names(crfModel@data@env$response)]])))^2
# set random seed
set.seed(seeds[iters])
# set random seed
set.seed(seeds[iters])
# conditional variable permutation importance
crfConPermImp <- permimp::permimp(crfModel, nperm=nperm, conditional=TRUE, threshold=permImpCondThres, progressBar=FALSE)
crfConPermImpVals <- crfConPermImpVals + crfConPermImp$values
crfConPermImpValsPerTree <- rbind(crfConPermImpValsPerTree, crfConPermImp$perTree)
}
# average metrics
crfOOBErrAll <- crfOOBErrAll/length(seeds)
crfOOBRMSE <- crfOOBRMSE/length(seeds)
crfOOBMAE <- crfOOBMAE/length(seeds)
crfOOBErrR2 <- crfOOBErrR2/length(seeds)
crfConPermImpVals <- data.frame(CondPermImp=sort(crfConPermImpVals/length(seeds), decreasing=TRUE))
#crfConPermImpValsMn <- sort(colMeans(crfConPermImpValsPerTree), decreasing=TRUE) # mean of conditional variable importance values per tree (same as crfConPermImpVals)
crfConPermImpValsQtl <- data.frame(apply(crfConPermImpValsPerTree, 2, quantile, probs=c(0.25, 0.50, 0.75)))
resultsOut <- list('OOB_RMSE'=crfOOBRMSE, 'OOB_MAE'=crfOOBMAE, 'Rsquared'=crfOOBErrR2, 'conditional_permimp'=crfConPermImpVals, 'conditional_permimp_perTree'=crfConPermImpValsPerTree, 'conditional_permimp_qtl'=crfConPermImpValsQtl)
return(resultsOut)
}
crfReg <- function(dataIn, iVars, dVar, seeds, ntree, mtry, permImpCondThres=0.95, minsplit=20, minbucket=7, nperm=1){
# initialise variables
crfOOBErrAll <- 0
crfOOBRMSE <- 0
crfOOBMAE <- 0
crfOOBErrR2 <- 0
crfMarPermImpVals <- 0
crfConPermImpVals <- 0
crfMarPermImpValsPerTree <- data.frame()
crfConPermImpValsPerTree <- data.frame()
# formula for regression
formVars <- reformulate(iVars, dVar)
for (iters in 1:length(seeds)){
# set random seed
set.seed(seeds[iters])
# train crf model
crfModel <- party::cforest(formVars, data=dataIn,
controls=party::cforest_unbiased(ntree=ntree,
mtry=mtry,
minsplit=minsplit,
minbucket=minbucket))
# conditional variable permutation importance
crfConPermImp <- permimp::permimp(crfModel, nperm=nperm, conditional=TRUE, threshold=permImpCondThres, progressBar=FALSE)
crfConPermImpVals <- crfConPermImp$values
if (iters == 1){
crfConPermImpVals1 <- data.frame(CondPermImp=sort(crfConPermImpVals, decreasing=TRUE))
crfConPermImpValsPerTree1 <- crfConPermImp$perTree
crfConPermImpValsQtl1 <- data.frame(apply(crfConPermImpValsPerTree1, 2, quantile, probs=c(0.25, 0.50, 0.75)))
# get OOB predictions
crfModelOOB <- predict(crfModel, OOB=TRUE, type='response')
# get OOB error
crfModelOOBErr <- as.numeric(as.matrix(as.numeric(as.matrix(crfModelOOB))
- as.numeric(as.matrix(crfModel@data@env$response[[names(crfModel@data@env$response)]]))))
# OOB RMSE, error quartiles and Rsquared
crfOOBRMSE <- sqrt(mean(crfModelOOBErr^2))
crfOOBMAE <- crfOOBMAE + mean(abs(crfModelOOBErr))
crfOOBErrR2 <- cor(as.numeric(as.matrix(crfModelOOB)),
as.numeric(as.matrix(crfModel@data@env$response[[names(crfModel@data@env$response)]])))^2
}
else{
crfConPermImpValsN <- data.frame(CondPermImp=sort(crfConPermImpVals, decreasing=TRUE))
nVarImpChecks <- min(max(sum(crfConPermImpVals1 >= crfConPermImpVals1$CondPermImp[1]*0.1),
sum(crfConPermImpValsN >= crfConPermImpValsN$CondPermImp[1]*0.1)), 4) # number of variable importance values with a value at least 10% of the highest importance
if (sum(rownames(crfConPermImpVals1)[1:nVarImpChecks] != rownames(crfConPermImpValsN)[1:nVarImpChecks]) > 0){
warning("Permutation importance rank order within 10% of max differs between seeds: increase number of trees ('ntree') or number of permutations ('nperm'), or subsample of features ('mtry')")
}
else{
resultsOut <- list('OOB_errors'=crfModelOOBErr, 'OOB_RMSE'=crfOOBRMSE, 'OOB_MAE'=crfOOBMAE, 'Rsquared'=crfOOBErrR2, 'conditional_permimp'=crfConPermImpVals1, 'conditional_permimp_perTree'=crfConPermImpValsPerTree1, 'conditional_permpimp_qtl'=crfConPermImpValsQtl1)
return(resultsOut)
}
}
}
}
mtryTune <- function(dataIn, iVars, dVar, seed, ntrees, minsplit=20, minbucket=7){
formVars <- reformulate(iVars, dVar)
# set mtry values and corresponding iVars/mtry ratios
iVars_mtrys <- c(10.5, 5.25, 3.5, 2.75, 2.25, 1.75, 1.5, 1.25)
mtrys <- as.integer(length(iVars)/iVars_mtrys)
iVars_mtrys <- iVars_mtrys[mtrys >= 2] # remove 0 or 1 values
mtrys <- mtrys[mtrys >= 2] # remove 0 or 1 values
# remove any duplicated values
iVars_mtrys <- iVars_mtrys[!(duplicated(mtrys) | duplicated(mtrys, fromLast = TRUE))]
mtrys <- mtrys[!(duplicated(mtrys) | duplicated(mtrys, fromLast = TRUE))]
# ensure mtrys is less than length(iVars)
iVars_mtrys <- iVars_mtrys[mtrys < length(iVars)]
mtrys <- mtrys[mtrys < length(iVars)]
resRMSEMap = matrix(data=NA, nrow=length(mtrys), ncol=length(ntrees))
resRsquaredMap = matrix(data=NA, nrow=length(mtrys), ncol=length(ntrees))
resMAEMap = matrix(data=NA, nrow=length(mtrys), ncol=length(ntrees))
for (ii in seq(1, length(ntrees))){
set.seed(seed)
ntree = ntrees[ii]
tuneMod <- caret::train(formVars,
data=dataIn,
method='cforest',
controls=party::cforest_unbiased(ntree=ntree,
minsplit=minsplit,
minbucket=minbucket),
tuneGrid=data.frame(.mtry=mtrys),
trControl = trainControl(method = "oob"))
tuneMod$results <- cbind(tuneMod$results, data.frame(iVars_mtrys=iVars_mtrys))
print(tuneMod$results)
print(tuneMod$bestTune)
resRMSEMap[, ii] = tuneMod$results$RMSE
resRsquaredMap[, ii] = tuneMod$results$Rsquared
resMAEMap[, ii] = tuneMod$results$MAE
}
# convert to data frames with mtry as row names and ntree as column names, and convert to long format using tidyverse
resdfRMSEMap <- as.data.frame(resRMSEMap)
rownames(resdfRMSEMap) <- mtrys
colnames(resdfRMSEMap) <- ntrees
resdfRsquaredMap <- as.data.frame(resRsquaredMap)
rownames(resdfRsquaredMap) <- mtrys
colnames(resdfRsquaredMap) <- ntrees
resdfMAEMap <- as.data.frame(resMAEMap)
rownames(resdfMAEMap) <- mtrys
colnames(resdfMAEMap) <- ntrees
# convert dataframes to long format using tidyverse
resdfRMSEMap <- resdfRMSEMap |>
rownames_to_column('mtry') |>
gather(key='ntree', value='RMSE', -mtry)
resdfRsquaredMap <- resdfRsquaredMap |>
rownames_to_column('mtry') |>
gather(key='ntree', value='Rsquared', -mtry)
resdfMAEMap <- resdfMAEMap |>
rownames_to_column('mtry') |>
gather(key='ntree', value='MAE', -mtry)
# ensure ntree and mtry columns are ordered factors
resdfRMSEMap$ntree <- factor(resdfRMSEMap$ntree, levels=as.character(ntrees))
resdfRMSEMap$mtry <- factor(resdfRMSEMap$mtry, levels=as.character(mtrys))
resdfRsquaredMap$ntree <- factor(resdfRsquaredMap$ntree, levels=as.character(ntrees))
resdfRsquaredMap$mtry <- factor(resdfRsquaredMap$mtry, levels=as.character(mtrys))
resdfMAEMap$ntree <- factor(resdfMAEMap$ntree, levels=as.character(ntrees))
resdfMAEMap$mtry <- factor(resdfMAEMap$mtry, levels=as.character(mtrys))
# plot heatmaps using ggplot, with extreme (min or max) value plotted as overlaid point using annotate and colourbar scale reversed
pHeatmapRMSE <- ggplot(resdfRMSEMap) +
geom_tile(aes(x=ntree, y=mtry, fill=RMSE)) +
scale_fill_viridis(option="viridis", direction=-1) +
geom_point(data=resdfRMSEMap[which(resdfRMSEMap$RMSE
== min(resdfRMSEMap$RMSE),
arr.ind = TRUE),],
aes(x=ntree, y=mtry), colour="red", size=2) +
guides(colour = guide_colourbar(reverse=TRUE)) +
labs(x="ntree", y="mtry", fill="RMSE") +
theme(text = element_text(family = "serif"))
pHeatmapRsquared <- ggplot(resdfRsquaredMap) +
geom_tile(aes(x=ntree, y=mtry, fill=Rsquared)) +
scale_fill_viridis(option="viridis", direction=1) +
geom_point(data=resdfRsquaredMap[which(resdfRsquaredMap$Rsquared
== max(resdfRsquaredMap$Rsquared),
arr.ind = TRUE),],
aes(x=ntree, y=mtry), colour="red", size=2) +
guides(colour = guide_colourbar(reverse=TRUE)) +
labs(x="ntree", y="mtry", fill="Rsquared") +
theme(text = element_text(family = "serif"))
pHeatmapMAE <- ggplot(resdfMAEMap) +
geom_tile(aes(x=ntree, y=mtry, fill=MAE)) +
scale_fill_viridis(option="viridis", direction=-1) +
geom_point(data=resdfMAEMap[which(resdfMAEMap$MAE
== min(resdfMAEMap$MAE),
arr.ind = TRUE),],
aes(x=ntree, y=mtry), colour="red", size=2) +
guides(colour = guide_colourbar(reverse=TRUE)) +
labs(x="ntree", y="mtry", fill="MAE") +
theme(text = element_text(family = "serif"))
p <- cowplot::plot_grid(pHeatmapRMSE,
pHeatmapRsquared,
pHeatmapMAE,
ncol=3, nrow=1)
return(p)
} # end of function
permImpCondThres <- 0.95
minsplit <- 20
minbucket <- 7
ntrees <- c(251, 501, 1001, 1501, 2501, 4001, 5501)
Initialise results output variables
resdAnnoyMnFit <- data.frame(RMSE = numeric(),
MAE = numeric(),
Rsquared = numeric())
resdAnnoyMnPermImp <- list()
Here, the ‘ambient only’ stimuli are removed, as the analysis cannot handle missing values for dB metrics.
stimDataANum <- stimDataANum[complete.cases(stimDataANum),]
stimDataANum$UASLAeq <- as.numeric(stimDataANum$UASLAeq)
stimDataANum$SNRlevel <- as.numeric(stimDataANum$SNRlevel)
iVars <- names(stimDataANum)[which(names(stimDataANum) == 'UASLAeq'):which(names(stimDataANum) == 'UASImpulsSHM05ExMaxLR')]
iVars <- iVars[! iVars %in% 'SNRlevel']
iVars <- c(iVars,
names(stimDataANum)[which(colnames(stimDataANum)=="LAeqLAF90diff"):
which(colnames(stimDataANum)=="dImpulsSHM05ExMaxLR")], "SNRlevel")
dVar <- "dAnnoyMean"
seeds <- c(14569, 98651, 54654498, 454948, 41321)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
if (saveplots){
ggsave(filename="PtAdAnnoyMnAllVarsHyperTune.svg", width=12, height=4, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnAllVarsHyperTune.svg")
ggsave(filename="PtAdAnnoyMnAllVarsHyperTune.pdf", width=12, height=4, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnAllVarsHyperTune.pdf")
}
Selected hyperparameters
ntree <- 251
mtry <- as.integer(length(iVars)/1.5)
Train preliminary model
nperm <- 10
resultsOutAbsDiffs <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutAbsDiffs$OOB_RMSE
[1] 0.434307
resultsOutAbsDiffs$OOB_MAE
[1] 0.3533671
resultsOutAbsDiffs$Rsquared
[1] 0.8738479
Train multiple seeds model
resultsOutAbsDiffs <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutAbsDiffs$OOB_RMSE
[1] 0.430126
resultsOutAbsDiffs$OOB_MAE
[1] 0.3537057
resultsOutAbsDiffs$Rsquared
[1] 0.875829
# store results
resdAnnoyMnFit['All vars', 'RMSE'] <- resultsOutAbsDiffs$OOB_RMSE
resdAnnoyMnFit['All vars', 'MAE'] <- resultsOutAbsDiffs$OOB_MAE
resdAnnoyMnFit['All vars', 'Rsquared'] <- resultsOutAbsDiffs$Rsquared
resdAnnoyMnPermImp$AllVars <- resultsOutAbsDiffs$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutAbsDiffs.conimp <- arrange(resultsOutAbsDiffs$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutAbsDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimp), levels=rownames(resultsOutAbsDiffs.conimp)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnAllVarsConPermimp.svg", width=8, height=20, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnAllVarsConPermimp.svg")
ggsave(filename="PtAdAnnoyMnAllVarsConPermimp.pdf", width=8, height=20, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnAllVarsConPermimp.pdf")
}
# Plot only positive values
resultsOutAbsDiffs.conimpPtv <- resultsOutAbsDiffs.conimp |>
rownames_to_column('Metric') |>
filter_if(is.numeric, all_vars(. > 0)) |>
column_to_rownames('Metric')
pBar <- ggplot(resultsOutAbsDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimpPtv), levels=rownames(resultsOutAbsDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnAllVarsConPermimpPtv.svg", width=8, height=8, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnAllVarsConPermimp.svg")
ggsave(filename="PtAdAnnoyMnAllVarsConPermimpPtv.pdf", width=8, height=8, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnAllVarsConPermimp.pdf")
}
iVars <- names(stimDataANum)[which(names(stimDataANum) == 'UASLAeq'):which(names(stimDataANum) == 'UASImpulsSHM05ExMaxLR')]
iVars <- iVars[! iVars %in% c('SNRlevel')]
dVar <- "dAnnoyMean"
seeds <- c(578312, 544, 84894, 54654, 153157)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
if (saveplots){
ggsave(filename="PtAdAnnoyMnAbsVarsHyperTune.svg", width=12, height=4, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnAbsVarsHyperTune.svg")
ggsave(filename="PtAdAnnoyMnAbsVarsHyperTune.pdf", width=12, height=4, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnAbsVarsHyperTune.pdf")
}
Selected hyperparameters
ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutAbs <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutAbs$OOB_RMSE
[1] 0.5920291
resultsOutAbs$OOB_MAE
[1] 0.4797379
resultsOutAbs$Rsquared
[1] 0.8046045
Train multiple seeds model
resultsOutAbs <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutAbs$OOB_RMSE
[1] 0.5898948
resultsOutAbs$OOB_MAE
[1] 0.4771605
resultsOutAbs$Rsquared
[1] 0.8064115
# store results
resdAnnoyMnFit['Abs vars', 'RMSE'] <- resultsOutAbs$OOB_RMSE
resdAnnoyMnFit['Abs vars', 'MAE'] <- resultsOutAbs$OOB_MAE
resdAnnoyMnFit['Abs vars', 'Rsquared'] <- resultsOutAbs$Rsquared
resdAnnoyMnPermImp$AbsVars <- resultsOutAbs$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutAbs.conimp <- arrange(resultsOutAbs$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutAbs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimp), levels=rownames(resultsOutAbs.conimp)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) +
coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnAbsVarsConPermimp.svg", width=8, height=11.5, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnAbsVarsConPermimp.svg")
ggsave(filename="PtAdAnnoyMnAbsVarsConPermimp.pdf", width=8, height=11.5, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnAbsVarsConPermimp.pdf")
}
# Plot values disregarding ambient LAeq
pBar <- ggplot(slice(resultsOutAbs.conimp, 1:n()-1)) + geom_col(aes(x=factor(rownames(slice(resultsOutAbs.conimp, 1:n()-1)), levels=rownames(slice(resultsOutAbs.conimp, 1:n()-1))), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnAbsVarsSkipAmbConPermimp.svg", width=8, height=11.5, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnAbsVarsSkipAmbConPermimp.svg")
ggsave(filename="PtAdAnnoyMnAbsVarsSkipAmbConPermimp.pdf", width=8, height=11.5, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnAbsVarsSkipAmbConPermimp.pdf")
}
# Plot only positive values disregarding ambient LAeq
resultsOutAbs.conimpPtv <- resultsOutAbs.conimp |>
rownames_to_column('Metric') |>
filter_if(is.numeric, all_vars(. > 0)) |>
column_to_rownames('Metric')
pBar <- ggplot(slice(resultsOutAbs.conimpPtv, 1:n()-1)) + geom_col(aes(x=factor(rownames(slice(resultsOutAbs.conimpPtv, 1:n()-1)), levels=rownames(slice(resultsOutAbs.conimpPtv, 1:n()-1))), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnAbsVarsSkipAmbConPermimpPtv.svg", width=8, height=7, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnAbsVarsSkipAmbConPermimpPtv.svg")
ggsave(filename="PtAdAnnoyMnAbsVarsSkipAmbConPermimpPtv.pdf", width=8, height=7, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnAbsVarsSkipAmbConPermimpPtv.pdf")
}
Selected metric
absVar <- "UASLoudECMAPowAvgBin"
Here, the ‘ambient only’ stimuli are replaced, as the SQM analysis can incorporate 0 values (whereas dB metrics cannot be handled in this way).
stimDataANum <- rbind(stimDataANum, stimDataA[stimDataA['Stimulus']
== "A1",
colnames(stimDataANum)],
stimDataA[stimDataA['Stimulus']
== "A2",
colnames(stimDataANum)])
stimDataANum <- arrange(stimDataANum, Stimulus)
iVars <- c(absVar, "AmbientLAeq", "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonLdECMAPowAvgBin", "UASTonLdECMA05ExBin", "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR")
dVar <- "dAnnoyMean"
seeds <- c(540, 104798, 456464, 87331, 94564)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 1001
mtry <- as.integer(length(iVars)/1.5)
Train preliminary model
# Tonality with tonal loudness
nperm <- 5
resultsOutTonal1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal1$OOB_RMSE
[1] 0.6390897
resultsOutTonal1$OOB_MAE
[1] 0.5098136
resultsOutTonal1$Rsquared
[1] 0.7782498
Train multiple seeds model
# Tonality with tonal loudness
resultsOutTonal1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal1$OOB_RMSE
[1] 0.6330849
resultsOutTonal1$OOB_MAE
[1] 0.5045541
resultsOutTonal1$Rsquared
[1] 0.7848499
# store results
resdAnnoyMnFit['Abs tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resdAnnoyMnFit['Abs tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resdAnnoyMnFit['Abs tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resdAnnoyMnPermImp$AbsTonal1 <- resultsOutTonal1$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Tonality inc. tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.2))
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnTonalLdConPermimp.svg", width=8, height=3.8, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnTonalLdConPermimp.svg")
ggsave(filename="PtAdAnnoyMnTonalLdConPermimp.pdf", width=8, height=3.8, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnTonalLdConPermimp.pdf")
}
Selected metric
tonLdVar <- "UASTonLdECMAPowAvgBin"
iVars <- c(absVar, "AmbientLAeq", "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR")
dVar <- "dAnnoyMean"
seeds <- c(156089, 5860, 10528, 89541, 4685146)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 1001
mtry <- as.integer(length(iVars)/1.5)
Train preliminary model
# Tonality
nperm <- 5
resultsOutTonal2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal2$OOB_RMSE
[1] 0.6030948
resultsOutTonal2$OOB_MAE
[1] 0.4825113
resultsOutTonal2$Rsquared
[1] 0.8186149
Train multiple seeds model
# Tonality
resultsOutTonal2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal2$OOB_RMSE
[1] 0.6032431
resultsOutTonal2$OOB_MAE
[1] 0.4815779
resultsOutTonal2$Rsquared
[1] 0.8154112
# store results
resdAnnoyMnFit['Abs tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resdAnnoyMnFit['Abs tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resdAnnoyMnFit['Abs tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resdAnnoyMnPermImp$AbsTonal2 <- resultsOutTonal2$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Tonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.2))
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnTonalConPermimp.svg", width=8, height=3.2, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnTonalConPermimp.svg")
ggsave(filename="PtAdAnnoyMnTonalConPermimp.pdf", width=8, height=3.2, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnTonalConPermimp.pdf")
}
Selected metric
tonalVar <- "UASTonalSHMInt05ExMaxLR"
# Fluctuation strength
iVars <- c(absVar, "AmbientLAeq", "UASFluctSHM10ExBin", "UASFluctSHM05ExBin", "UASFluctFZ10ExMaxLR", "UASFluctFZ05ExMaxLR", "UASFluctOV10ExMaxLR", "UASFluctOV05ExMaxLR")
dVar <- "dAnnoyMean"
seeds <- c(25107, 546098, 195, 5937, 102658)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutFluct <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutFluct$OOB_RMSE
[1] 0.5843817
resultsOutFluct$OOB_MAE
[1] 0.4552687
resultsOutFluct$Rsquared
[1] 0.8202794
Train multiple seeds model
resultsOutFluct <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutFluct$OOB_RMSE
[1] 0.5847776
resultsOutFluct$OOB_MAE
[1] 0.456084
resultsOutFluct$Rsquared
[1] 0.8194579
# store results
resdAnnoyMnFit['Abs fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resdAnnoyMnFit['Abs fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resdAnnoyMnFit['Abs fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resdAnnoyMnPermImp$AbsFluct <- resultsOutFluct$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Fluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnFluctConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnFluctConPermimp.svg")
ggsave(filename="PtAdAnnoyMnFluctConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnFluctConPermimp.pdf")
}
Selected metric
fluctVar <- "UASFluctOV10ExMaxLR"
# Roughness
iVars <- c(absVar, "AmbientLAeq", "UASRoughECMA10ExBin", "UASRoughECMA05ExBin", "UASRoughFZ10ExMaxLR", "UASRoughFZ05ExMaxLR", "UASRoughDW10ExMaxLR", "UASRoughDW05ExMaxLR")
dVar <- "dAnnoyMean"
seeds <- c(4701, 52187, 16589, 65217, 16893)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutRough <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutRough$OOB_RMSE
[1] 0.5978974
resultsOutRough$OOB_MAE
[1] 0.465474
resultsOutRough$Rsquared
[1] 0.8040179
Train multiple seeds model
resultsOutRough <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutRough$OOB_RMSE
[1] 0.5997308
resultsOutRough$OOB_MAE
[1] 0.4671883
resultsOutRough$Rsquared
[1] 0.8028445
# store results
resdAnnoyMnFit['Abs rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resdAnnoyMnFit['Abs rough', 'MAE'] <- resultsOutRough$OOB_MAE
resdAnnoyMnFit['Abs rough', 'Rsquared'] <- resultsOutRough$Rsquared
resdAnnoyMnPermImp$AbsRough <- resultsOutRough$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Roughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnRoughConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnRoughConPermimp.svg")
ggsave(filename="PtAdAnnoyMnRoughConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnRoughConPermimp.pdf")
}
Selected metric
roughVar <- "UASRoughFZ05ExMaxLR"
# Impulsiveness
iVars <- c(absVar, "AmbientLAeq", "UASImpulsSHMAvgMaxLR", "UASImpulsSHM05ExMaxLR", "UASImpulsSHMPowAvgMaxLR")
dVar <- "dAnnoyMean"
seeds <- c(8495, 59867, 5416, 9843, 86)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 251
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutImpuls <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutImpuls$OOB_RMSE
[1] 0.6105395
resultsOutImpuls$OOB_MAE
[1] 0.4733032
resultsOutImpuls$Rsquared
[1] 0.7863911
Train multiple seeds model
resultsOutImpuls <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutImpuls$OOB_RMSE
[1] 0.603079
resultsOutImpuls$OOB_MAE
[1] 0.4665693
resultsOutImpuls$Rsquared
[1] 0.7960339
# store results
resdAnnoyMnFit['Abs impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resdAnnoyMnFit['Abs impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resdAnnoyMnFit['Abs impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resdAnnoyMnPermImp$AbsImpuls <- resultsOutImpuls$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("Impulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnImpulsConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnImpulsConPermimp.svg")
ggsave(filename="PtAdAnnoyMnImpulsConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnImpulsConPermimp.pdf")
}
Selected metric
impulsVar <- "UASImpulsSHMAvgMaxLR"
Now the highest importance SQMs are ranked against each other, controlling for UAS loudness and ambient LAeq.
iVars <- c(absVar, "AmbientLAeq", sharpVar, tonLdVar, fluctVar, roughVar, impulsVar)
dVar <- "dAnnoyMean"
seeds <- c(70498, 4, 14986, 453, 864)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 1001
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutSQMs1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs1$OOB_RMSE
[1] 0.6099365
resultsOutSQMs1$OOB_MAE
[1] 0.4869588
resultsOutSQMs1$Rsquared
[1] 0.8025851
Train multiple seeds model
resultsOutSQMs1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs1$OOB_RMSE
[1] 0.60767
resultsOutSQMs1$OOB_MAE
[1] 0.4855243
resultsOutSQMs1$Rsquared
[1] 0.805056
# store results
resdAnnoyMnFit['Abs SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resdAnnoyMnFit['Abs SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resdAnnoyMnFit['Abs SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resdAnnoyMnPermImp$AbsSQMs1 <- resultsOutSQMs1$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1))
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnAbsSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnAbsSQMsTonLdConPermimp.svg")
ggsave(filename="PtAdAnnoyMnAbsSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnAbsSQMsTonLdConPermimp.pdf")
}
iVars <- c(absVar, "AmbientLAeq", sharpVar, tonalVar, fluctVar, roughVar, impulsVar)
dVar <- "dAnnoyMean"
seeds <- c(546, 57203, 270835, 60592, 8094)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 4001
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutSQMs2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs2$OOB_RMSE
[1] 0.5870765
resultsOutSQMs2$OOB_MAE
[1] 0.4689816
resultsOutSQMs2$Rsquared
[1] 0.8254522
Train multiple seeds model
resultsOutSQMs2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs2$OOB_RMSE
[1] 0.5880773
resultsOutSQMs2$OOB_MAE
[1] 0.469515
resultsOutSQMs2$Rsquared
[1] 0.8255209
# store results
resdAnnoyMnFit['Abs SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resdAnnoyMnFit['Abs SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_MAE
resdAnnoyMnFit['Abs SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resdAnnoyMnPermImp$AbsSQMs2 <- resultsOutSQMs2$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1))
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnAbsSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnAbsSQMsNoTonLdConPermimp.svg")
ggsave(filename="PtAdAnnoyMnAbsSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnAbsSQMsNoTonLdConPermimp.pdf")
}
Next, the difference metrics are analysed,
Here, the ‘ambient only’ stimuli are removed, as the analysis cannot handle missing values for dB metrics.
stimDataANum <- stimDataANum[complete.cases(stimDataANum),]
stimDataANum$UASLAeq <- as.numeric(stimDataANum$UASLAeq)
stimDataANum$SNRlevel <- as.numeric(stimDataANum$SNRlevel)
iVars <- c(names(stimDataANum)[which(colnames(stimDataANum)=="LAeqLAF90diff"):
which(colnames(stimDataANum)=="dImpulsSHM05ExMaxLR")], "SNRlevel")
dVar <- "dAnnoyMean"
seeds <- c(568392, 498, 4089, 78132, 741809)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 251
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutDiffs <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutDiffs$OOB_RMSE
[1] 0.4308342
resultsOutDiffs$OOB_MAE
[1] 0.3567504
resultsOutDiffs$Rsquared
[1] 0.8740203
Train multiple seeds model
resultsOutDiffs <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutDiffs$OOB_RMSE
[1] 0.4314902
resultsOutDiffs$OOB_MAE
[1] 0.3566699
resultsOutDiffs$Rsquared
[1] 0.8741682
# store results
resdAnnoyMnFit['Diff vars', 'RMSE'] <- resultsOutDiffs$OOB_RMSE
resdAnnoyMnFit['Diff vars', 'MAE'] <- resultsOutDiffs$OOB_MAE
resdAnnoyMnFit['Diff vars', 'Rsquared'] <- resultsOutDiffs$Rsquared
resdAnnoyMnPermImp$DiffVars <- resultsOutDiffs$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutDiffs.conimp <- arrange(resultsOutDiffs$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimp), levels=rownames(resultsOutDiffs.conimp)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnDiffVarsConPermimp.svg", width=8, height=10, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnDiffVarsConPermimp.svg")
ggsave(filename="PtAdAnnoyMnDiffVarsConPermimp.pdf", width=8, height=10, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnDiffVarsConPermimp.pdf")
}
# Plot only positive values
resultsOutDiffs.conimpPtv <- resultsOutDiffs.conimp |>
rownames_to_column('Metric') |>
filter_if(is.numeric, all_vars(. > 0)) |>
column_to_rownames('Metric')
pBar <- ggplot(resultsOutDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimpPtv), levels=rownames(resultsOutDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnDiffVarsConPermimpPtv.svg", width=8, height=8, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnDiffVarsConPermimpPtv.svg")
ggsave(filename="PtAdAnnoyMnDiffVarsConPermimpPtv.pdf", width=8, height=8, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnDiffVarsConPermimpPtv.pdf")
}
Selected metric
diffVar <- "EPNLLAF50diff"
The ambient-only stimuli are NOT replaced here, as the difference analysis includes dB metrics.
iVars <- c(diffVar, "dSharpAurISO3PowAvgMaxLR", "dSharpAurISO305ExMaxLR", "dSharpAurSHMPowAvgMaxLR", "dSharpAurSHM05ExMaxLR")
dVar <- "dAnnoyMean"
seeds <- c(84194, 905, 64815, 928054, 625091)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutSharp <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSharp$OOB_RMSE
[1] 0.4523884
resultsOutSharp$OOB_MAE
[1] 0.3595752
resultsOutSharp$Rsquared
[1] 0.8627182
Train multiple seeds model
resultsOutSharp <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSharp$OOB_RMSE
[1] 0.4536154
resultsOutSharp$OOB_MAE
[1] 0.3606428
resultsOutSharp$Rsquared
[1] 0.8618904
# store results
resdAnnoyMnFit['Diff sharp', 'RMSE'] <- resultsOutSharp$OOB_RMSE
resdAnnoyMnFit['Diff sharp', 'MAE'] <- resultsOutSharp$OOB_MAE
resdAnnoyMnFit['Diff sharp', 'Rsquared'] <- resultsOutSharp$Rsquared
resdAnnoyMnPermImp$DiffSharp <- resultsOutSharp$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSharp.conimp <- arrange(resultsOutSharp$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSharp.conimp) + geom_col(aes(x=factor(rownames(resultsOutSharp.conimp), levels=rownames(resultsOutSharp.conimp)), y=CondPermImp), fill=mycolours[2], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dSharpness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMndSharpConPermimp.svg", width=8, height=2.6, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMndSharpConPermimp.svg")
ggsave(filename="PtAdAnnoyMndSharpConPermimp.pdf", width=8, height=2.6, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMndSharpConPermimp.pdf")
}
Selected metric
dSharpVar <- "dSharpAurISO3PowAvgMaxLR"
iVars <- c(diffVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR", "dTonLdECMAPowAvgBin", "dTonLdECMA05ExBin")
dVar <- "dAnnoyMean"
seeds <- c(561684, 104798, 1536, 48, 48561)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
# Tonality with tonal loudness
nperm <- 5
resultsOutTonal1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal1$OOB_RMSE
[1] 0.4541254
resultsOutTonal1$OOB_MAE
[1] 0.3676218
resultsOutTonal1$Rsquared
[1] 0.8616272
Train multiple seeds model
# Tonality with tonal loudness
resultsOutTonal1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal1$OOB_RMSE
[1] 0.4527969
resultsOutTonal1$OOB_MAE
[1] 0.3685465
resultsOutTonal1$Rsquared
[1] 0.8625633
# store results
resdAnnoyMnFit['Diff tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resdAnnoyMnFit['Diff tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resdAnnoyMnFit['Diff tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resdAnnoyMnPermImp$DiffTonal1 <- resultsOutTonal1$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dTonality inc. dtonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.6))
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMndTonalLdConPermimp.svg", width=8, height=2.6, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMndTonalLdConPermimp.svg")
ggsave(filename="PtAdAnnoyMndTonalLdConPermimp.pdf", width=8, height=2.6, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMndTonalLdConPermimp.pdf")
}
Selected metric
dTonLdVar <- "dTonLdECMAPowAvgBin"
iVars <- c(diffVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR")
dVar <- "dAnnoyMean"
seeds <- c(410865, 2954, 70812, 203, 7984)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 4001
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
# Tonality
nperm <- 5
resultsOutTonal2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal2$OOB_RMSE
[1] 0.462319
resultsOutTonal2$OOB_MAE
[1] 0.3667896
resultsOutTonal2$Rsquared
[1] 0.8545831
Train multiple seeds model
# Tonality
resultsOutTonal2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal2$OOB_RMSE
[1] 0.4619735
resultsOutTonal2$OOB_MAE
[1] 0.3664541
resultsOutTonal2$Rsquared
[1] 0.8548939
# store results
resdAnnoyMnFit['Diff tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resdAnnoyMnFit['Diff tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resdAnnoyMnFit['Diff tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resdAnnoyMnPermImp$DiffTonal2 <- resultsOutTonal2$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dTonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.6))
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMndTonalConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMndTonalConPermimp.svg")
ggsave(filename="PtAdAnnoyMndTonalConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMndTonalConPermimp.pdf")
}
Selectec metric
dTonalVar <- "dTonalSHMIntAvgMaxLR"
# Fluctuation strength
iVars <- c(diffVar, "dFluctSHM10ExBin", "dFluctSHM05ExBin", "dFluctOV10ExMaxLR", "dFluctOV05ExMaxLR")
dVar <- "dAnnoyMean"
seeds <- c(418657, 84, 1630, 18659, 3687)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 1001
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutFluct <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutFluct$OOB_RMSE
[1] 0.4693879
resultsOutFluct$OOB_MAE
[1] 0.3674453
resultsOutFluct$Rsquared
[1] 0.8512241
Train multiple seeds model
resultsOutFluct <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutFluct$OOB_RMSE
[1] 0.4696581
resultsOutFluct$OOB_MAE
[1] 0.368869
resultsOutFluct$Rsquared
[1] 0.8511147
# store results
resdAnnoyMnFit['Diff fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resdAnnoyMnFit['Diff fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resdAnnoyMnFit['Diff fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resdAnnoyMnPermImp$DiffFluct <- resultsOutFluct$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dFluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMndFluctConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnFluctConPermimp.svg")
ggsave(filename="PtAdAnnoyMndFluctConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnFluctConPermimp.pdf")
}
Selected metric
dFluctVar <- "dFluctOV10ExMaxLR"
# Roughness
iVars <- c(diffVar, "dRoughECMA10ExBin", "dRoughECMA05ExBin", "dRoughFZ10ExMaxLR", "dRoughFZ05ExMaxLR")
dVar <- "dAnnoyMean"
seeds <- c(69851, 85109, 410986, 1563, 896)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutRough <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutRough$OOB_RMSE
[1] 0.4931997
resultsOutRough$OOB_MAE
[1] 0.372368
resultsOutRough$Rsquared
[1] 0.839582
Train multiple seeds model
resultsOutRough <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutRough$OOB_RMSE
[1] 0.4955747
resultsOutRough$OOB_MAE
[1] 0.3727289
resultsOutRough$Rsquared
[1] 0.8381762
# store results
resdAnnoyMnFit['Diff rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resdAnnoyMnFit['Diff rough', 'MAE'] <- resultsOutRough$OOB_MAE
resdAnnoyMnFit['Diff rough', 'Rsquared'] <- resultsOutRough$Rsquared
resdAnnoyMnPermImp$DiffRough <- resultsOutRough$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dRoughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMndRoughConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMndRoughConPermimp.svg")
ggsave(filename="PtAdAnnoyMndRoughConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMndRoughConPermimp.pdf")
}
Selected metric
dRoughVar <- "dRoughECMA05ExBin"
# Impulsiveness
iVars <- c(diffVar, "dImpulsSHMAvgMaxLR", "dImpulsSHM05ExMaxLR", "dImpulsSHMPowAvgMaxLR")
dVar <- "dAnnoyMean"
seeds <- c(418659, 7805, 38475, 65834, 1653)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutImpuls <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutImpuls$OOB_RMSE
[1] 0.5242683
resultsOutImpuls$OOB_MAE
[1] 0.4079044
resultsOutImpuls$Rsquared
[1] 0.8135546
Train multiple seeds model
resultsOutImpuls <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutImpuls$OOB_RMSE
[1] 0.5207893
resultsOutImpuls$OOB_MAE
[1] 0.4047467
resultsOutImpuls$Rsquared
[1] 0.8173652
# store results
resdAnnoyMnFit['Diff impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resdAnnoyMnFit['Diff impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resdAnnoyMnFit['Diff impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resdAnnoyMnPermImp$DiffImpuls <- resultsOutImpuls$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + ggtitle("dImpulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMndImpulsConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMndImpulsConPermimp.svg")
ggsave(filename="PtAdAnnoyMndImpulsConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMndImpulsConPermimp.pdf")
}
Selected metric
dImpulsVar <- "dImpulsSHM05ExMaxLR"
Now the highest importance dSQMs are ranked against each other, controlling for loudness difference.
iVars <- c(diffVar, dSharpVar, dTonLdVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "dAnnoyMean"
seeds <- c(98465, 54163, 6541, 36485, 849675)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 251
mtry <- as.integer(length(iVars)/1.75)
Train preliminary model
nperm <- 5
resultsOutSQMs1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs1$OOB_RMSE
[1] 0.4579162
resultsOutSQMs1$OOB_MAE
[1] 0.3597023
resultsOutSQMs1$Rsquared
[1] 0.8618785
Train multiple seeds model
resultsOutSQMs1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs1$OOB_RMSE
[1] 0.460312
resultsOutSQMs1$OOB_MAE
[1] 0.363146
resultsOutSQMs1$Rsquared
[1] 0.8597486
# store results
resdAnnoyMnFit['Diff SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resdAnnoyMnFit['Diff SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resdAnnoyMnFit['Diff SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resdAnnoyMnPermImp$DiffSQMs1 <- resultsOutSQMs1$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.2))
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnDiffSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnDiffSQMsTonLdConPermimp.svg")
ggsave(filename="PtAdAnnoyMnDiffSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnDiffSQMsTonLdConPermimp.pdf")
}
iVars <- c(diffVar, dSharpVar, dTonalVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "dAnnoyMean"
seeds <- c(49865, 7852, 845961, 410583, 36748)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 4001
mtry <- as.integer(length(iVars)/1.75)
Train preliminary model
nperm <- 5
resultsOutSQMs2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs2$OOB_RMSE
[1] 0.458182
resultsOutSQMs2$OOB_MAE
[1] 0.360887
resultsOutSQMs2$Rsquared
[1] 0.862422
Train multiple seeds model
resultsOutSQMs2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs2$OOB_RMSE
[1] 0.4569997
resultsOutSQMs2$OOB_MAE
[1] 0.3596057
resultsOutSQMs2$Rsquared
[1] 0.863165
# store results
resdAnnoyMnFit['Diff SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resdAnnoyMnFit['Diff SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_MAE
resdAnnoyMnFit['Diff SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resdAnnoyMnPermImp$DiffSQMs2 <- resultsOutSQMs2$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (mean change in annoyance)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1.2))
pBar
if (saveplots){
ggsave(filename="PtAdAnnoyMnDiffSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
unlink("PtAdAnnoyMnDiffSQMsNoTonLdConPermimp.svg")
ggsave(filename="PtAdAnnoyMnDiffSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
unlink("PtAdAnnoyMnDiffSQMsNoTonLdConPermimp.pdf")
}
if (savedata){
utils::write.csv(resdAnnoyMnFit, paste(outDataPath, "\\ptACRFdAnnoyMnOOBFit.csv", sep=""))
ii <- 0
temp = list()
for (res in resdAnnoyMnPermImp){
ii <- ii + 1
temp[[ii]] <- as.data.frame(resdAnnoyMnPermImp[ii])
names(temp[[ii]]) <- names(resdAnnoyMnPermImp[ii])
}
openxlsx::write.xlsx(temp, paste(outDataPath, "\\ptACRFdAnnoyMnConPermimp.xlsx",
sep=""),
rowNames=TRUE)
}
Initialise results output variables
resdHiAnnoyFit <- data.frame(RMSE = numeric(),
MAE = numeric(),
Rsquared = numeric())
resdHiAnnoyPermImp <- list()
Here, the ‘ambient only’ stimuli are removed, as the analysis cannot handle missing values for dB metrics.
stimDataANum <- stimDataANum[complete.cases(stimDataANum),]
stimDataANum$UASLAeq <- as.numeric(stimDataANum$UASLAeq)
stimDataANum$SNRlevel <- as.numeric(stimDataANum$SNRlevel)
iVars <- names(stimDataANum)[which(names(stimDataANum) == 'UASLAeq'):which(names(stimDataANum) == 'UASImpulsSHM05ExMaxLR')]
iVars <- iVars[! iVars %in% 'SNRlevel']
iVars <- c(iVars,
names(stimDataANum)[which(colnames(stimDataANum)=="LAeqLAF90diff"):
which(colnames(stimDataANum)=="dImpulsSHM05ExMaxLR")], "SNRlevel")
dVar <- "dHighAnnoyPc"
seeds <- c(2, 312, 1897, 465978, 821659)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 501
mtry <- as.integer(length(iVars)/1.75)
Train preliminary model
nperm <- 5
resultsOutAbsDiffs <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutAbsDiffs$OOB_RMSE
[1] 5.468804
resultsOutAbsDiffs$OOB_MAE
[1] 4.327614
resultsOutAbsDiffs$Rsquared
[1] 0.6347201
Train multiple seeds model
resultsOutAbsDiffs <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutAbsDiffs$OOB_RMSE
[1] 5.504329
resultsOutAbsDiffs$OOB_MAE
[1] 4.353992
resultsOutAbsDiffs$Rsquared
[1] 0.6301227
# store results
resdHiAnnoyFit['All vars', 'RMSE'] <- resultsOutAbsDiffs$OOB_RMSE
resdHiAnnoyFit['All vars', 'MAE'] <- resultsOutAbsDiffs$OOB_MAE
resdHiAnnoyFit['All vars', 'Rsquared'] <- resultsOutAbsDiffs$Rsquared
resdHiAnnoyPermImp$AllVars <- resultsOutAbsDiffs$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutAbsDiffs.conimp <- arrange(resultsOutAbsDiffs$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutAbsDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimp), levels=rownames(resultsOutAbsDiffs.conimp)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyAllVarsConPermimp.svg", width=8, height=20, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyAllVarsConPermimp.svg")
ggsave(filename="PtAdHiAnnoyAllVarsConPermimp.pdf", width=8, height=20, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyAllVarsConPermimp.pdf")
}
# Plot only positive values
resultsOutAbsDiffs.conimpPtv <- resultsOutAbsDiffs.conimp |>
rownames_to_column('Metric') |>
filter_if(is.numeric, all_vars(. > 0)) |>
column_to_rownames('Metric')
pBar <- ggplot(resultsOutAbsDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimpPtv), levels=rownames(resultsOutAbsDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyAllVarsConPermimpPtv.svg", width=8, height=12, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyAllVarsConPermimp.svg")
ggsave(filename="PtAdHiAnnoyAllVarsConPermimpPtv.pdf", width=8, height=12, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyAllVarsConPermimp.pdf")
}
iVars <- names(stimDataANum)[which(names(stimDataANum) == 'UASLAeq'):which(names(stimDataANum) == 'UASImpulsSHM05ExMaxLR')]
iVars <- iVars[! iVars %in% c('SNRlevel')]
dVar <- "dHighAnnoyPc"
seeds <- c(578312, 544, 84894, 54654, 153157)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
if (saveplots){
ggsave(filename="PtAdHiAnnoyAbsVarsHyperTune.svg", width=12, height=4, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyAbsVarsHyperTune.svg")
ggsave(filename="PtAdHiAnnoyAbsVarsHyperTune.pdf", width=12, height=4, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyAbsVarsHyperTune.pdf")
}
Selected hyperparameters
ntree <- 501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 20
resultsOutAbs <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutAbs$OOB_RMSE
[1] 5.451148
resultsOutAbs$OOB_MAE
[1] 4.24412
resultsOutAbs$Rsquared
[1] 0.6356777
Train multiple seeds model
resultsOutAbs <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutAbs$OOB_RMSE
[1] 5.450724
resultsOutAbs$OOB_MAE
[1] 4.258414
resultsOutAbs$Rsquared
[1] 0.6356609
# store results
resdHiAnnoyFit['Abs vars', 'RMSE'] <- resultsOutAbs$OOB_RMSE
resdHiAnnoyFit['Abs vars', 'MAE'] <- resultsOutAbs$OOB_MAE
resdHiAnnoyFit['Abs vars', 'Rsquared'] <- resultsOutAbs$Rsquared
resdHiAnnoyPermImp$AbsVars <- resultsOutAbs$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutAbs.conimp <- arrange(resultsOutAbs$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutAbs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimp), levels=rownames(resultsOutAbs.conimp)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) +
coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyAbsVarsConPermimp.svg", width=8, height=11.5, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyAbsVarsConPermimp.svg")
ggsave(filename="PtAdHiAnnoyAbsVarsConPermimp.pdf", width=8, height=11.5, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyAbsVarsConPermimp.pdf")
}
# Plot only positive values
resultsOutAbs.conimpPtv <- resultsOutAbs.conimp |>
rownames_to_column('Metric') |>
filter_if(is.numeric, all_vars(. > 0)) |>
column_to_rownames('Metric')
pBar <- ggplot(resultsOutAbs.conimpPtv,) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimpPtv), levels=rownames(resultsOutAbs.conimpPtv)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyAbsVarsConPermimpPtv.svg", width=8, height=7, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyAbsVarsConPermimpPtv.svg")
ggsave(filename="PtAdHiAnnoyAbsVarsConPermimpPtv.pdf", width=8, height=7, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyAbsVarsConPermimpPtv.pdf")
}
Selected metric
absVar <- "UASLoudISO3PowAvgBin"
Here, the ‘ambient only’ stimuli are replaced, as the SQM analysis can incorporate 0 values (whereas dB metrics cannot be handled in this way).
stimDataANum <- rbind(stimDataANum, stimDataA[stimDataA['Stimulus']
== "A1",
colnames(stimDataANum)],
stimDataA[stimDataA['Stimulus']
== "A2",
colnames(stimDataANum)])
stimDataANum <- arrange(stimDataANum, Stimulus)
iVars <- c(absVar, "AmbientLAeq", "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonLdECMAPowAvgBin", "UASTonLdECMA05ExBin", "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR")
dVar <- "dHighAnnoyPc"
seeds <- c(540, 104798, 456464, 87331, 94564)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 251
mtry <- as.integer(length(iVars)/1.5)
Train preliminary model
# Tonality with tonal loudness
nperm <- 5
resultsOutTonal1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal1$OOB_RMSE
[1] 5.49583
resultsOutTonal1$OOB_MAE
[1] 4.367014
resultsOutTonal1$Rsquared
[1] 0.6298765
Train multiple seeds model
# Tonality with tonal loudness
resultsOutTonal1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal1$OOB_RMSE
[1] 5.456376
resultsOutTonal1$OOB_MAE
[1] 4.328807
resultsOutTonal1$Rsquared
[1] 0.6351314
# store results
resdHiAnnoyFit['Abs tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resdHiAnnoyFit['Abs tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resdHiAnnoyFit['Abs tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resdHiAnnoyPermImp$AbsTonal1 <- resultsOutTonal1$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Tonality inc. tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 70))
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyTonalLdConPermimp.svg", width=8, height=3.8, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyTonalLdConPermimp.svg")
ggsave(filename="PtAdHiAnnoyTonalLdConPermimp.pdf", width=8, height=3.8, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyTonalLdConPermimp.pdf")
}
Selected metric
tonLdVar <- "UASTonLdECMAPowAvgBin"
iVars <- c(absVar, "AmbientLAeq", "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR")
dVar <- "dHighAnnoyPc"
seeds <- c(156089, 5860, 10528, 89541, 4685146)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
# Tonality
nperm <- 5
resultsOutTonal2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm,
minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal2$OOB_RMSE
[1] 5.450319
resultsOutTonal2$OOB_MAE
[1] 4.329233
resultsOutTonal2$Rsquared
[1] 0.6372464
Train multiple seeds model
# Tonality
resultsOutTonal2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal2$OOB_RMSE
[1] 5.438749
resultsOutTonal2$OOB_MAE
[1] 4.317752
resultsOutTonal2$Rsquared
[1] 0.6388582
# store results
resdHiAnnoyFit['Abs tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resdHiAnnoyFit['Abs tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resdHiAnnoyFit['Abs tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resdHiAnnoyPermImp$AbsTonal2 <- resultsOutTonal2$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Tonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 70))
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyTonalConPermimp.svg", width=8, height=3.2, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyTonalConPermimp.svg")
ggsave(filename="PtAdHiAnnoyTonalConPermimp.pdf", width=8, height=3.2, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyTonalConPermimp.pdf")
}
Selected metric
tonalVar <- "UASTonalSHMInt05ExMaxLR"
# Fluctuation strength
iVars <- c(absVar, "AmbientLAeq", "UASFluctSHM10ExBin", "UASFluctSHM05ExBin", "UASFluctFZ10ExMaxLR", "UASFluctFZ05ExMaxLR", "UASFluctOV10ExMaxLR", "UASFluctOV05ExMaxLR")
dVar <- "dHighAnnoyPc"
seeds <- c(25107, 546098, 195, 5937, 102658)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutFluct <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2],
ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres,
nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutFluct$OOB_RMSE
[1] 5.663893
resultsOutFluct$OOB_MAE
[1] 4.453616
resultsOutFluct$Rsquared
[1] 0.6119349
Train multiple seeds model
resultsOutFluct <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutFluct$OOB_RMSE
[1] 5.631703
resultsOutFluct$OOB_MAE
[1] 4.416416
resultsOutFluct$Rsquared
[1] 0.6156464
# store results
resdHiAnnoyFit['Abs fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resdHiAnnoyFit['Abs fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resdHiAnnoyFit['Abs fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resdHiAnnoyPermImp$AbsFluct <- resultsOutFluct$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Fluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyFluctConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyFluctConPermimp.svg")
ggsave(filename="PtAdHiAnnoyFluctConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyFluctConPermimp.pdf")
}
Selected metric
fluctVar <- "UASFluctOV10ExMaxLR"
# Roughness
iVars <- c(absVar, "AmbientLAeq", "UASRoughECMA10ExBin", "UASRoughECMA05ExBin", "UASRoughFZ10ExMaxLR", "UASRoughFZ05ExMaxLR", "UASRoughDW10ExMaxLR", "UASRoughDW05ExMaxLR")
dVar <- "dHighAnnoyPc"
seeds <- c(4701, 52187, 16589, 65217, 16893)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutRough <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutRough$OOB_RMSE
[1] 5.582484
resultsOutRough$OOB_MAE
[1] 4.337246
resultsOutRough$Rsquared
[1] 0.6209775
Train multiple seeds model
resultsOutRough <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutRough$OOB_RMSE
[1] 5.583989
resultsOutRough$OOB_MAE
[1] 4.331894
resultsOutRough$Rsquared
[1] 0.621023
# store results
resdHiAnnoyFit['Abs rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resdHiAnnoyFit['Abs rough', 'MAE'] <- resultsOutRough$OOB_MAE
resdHiAnnoyFit['Abs rough', 'Rsquared'] <- resultsOutRough$Rsquared
resdHiAnnoyPermImp$AbsRough <- resultsOutRough$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Roughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyRoughConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyRoughConPermimp.svg")
ggsave(filename="PtAdHiAnnoyRoughConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyRoughConPermimp.pdf")
}
Selected metric
roughVar <- "UASRoughFZ05ExMaxLR"
# Impulsiveness
iVars <- c(absVar, "AmbientLAeq", "UASImpulsSHMAvgMaxLR", "UASImpulsSHM05ExMaxLR", "UASImpulsSHMPowAvgMaxLR")
dVar <- "dHighAnnoyPc"
seeds <- c(8495, 59867, 5416, 9843, 86)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutImpuls <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutImpuls$OOB_RMSE
[1] 5.494835
resultsOutImpuls$OOB_MAE
[1] 4.311435
resultsOutImpuls$Rsquared
[1] 0.6334893
Train multiple seeds model
resultsOutImpuls <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutImpuls$OOB_RMSE
[1] 5.482274
resultsOutImpuls$OOB_MAE
[1] 4.301147
resultsOutImpuls$Rsquared
[1] 0.6348031
# store results
resdHiAnnoyFit['Abs impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resdHiAnnoyFit['Abs impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resdHiAnnoyFit['Abs impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resdHiAnnoyPermImp$AbsImpuls <- resultsOutImpuls$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("Impulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyImpulsConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyImpulsConPermimp.svg")
ggsave(filename="PtAdHiAnnoyImpulsConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyImpulsConPermimp.pdf")
}
Selected metric
impulsVar <- "UASImpulsSHMAvgMaxLR"
Now the highest importance SQMs are ranked against each other, controlling for UAS loudness and ambient LAeq.
iVars <- c(absVar, "AmbientLAeq", sharpVar, tonLdVar, fluctVar, roughVar, impulsVar)
dVar <- "dHighAnnoyPc"
seeds <- c(70498, 4, 14986, 453, 864)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutSQMs1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs1$OOB_RMSE
[1] 5.422008
resultsOutSQMs1$OOB_MAE
[1] 4.293746
resultsOutSQMs1$Rsquared
[1] 0.6411011
Train multiple seeds model
resultsOutSQMs1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs1$OOB_RMSE
[1] 5.424953
resultsOutSQMs1$OOB_MAE
[1] 4.297558
resultsOutSQMs1$Rsquared
[1] 0.6407472
# store results
resdHiAnnoyFit['Abs SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resdHiAnnoyFit['Abs SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resdHiAnnoyFit['Abs SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resdHiAnnoyPermImp$AbsSQMs1 <- resultsOutSQMs1$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 50))
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyAbsSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyAbsSQMsTonLdConPermimp.svg")
ggsave(filename="PtAdHiAnnoyAbsSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyAbsSQMsTonLdConPermimp.pdf")
}
iVars <- c(absVar, "AmbientLAeq", sharpVar, tonalVar, fluctVar, roughVar, impulsVar)
dVar <- "dHighAnnoyPc"
seeds <- c(546, 57203, 270835, 60592, 8094)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutSQMs2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs2$OOB_RMSE
[1] 5.541448
resultsOutSQMs2$OOB_MAE
[1] 4.390907
resultsOutSQMs2$Rsquared
[1] 0.6273827
Train multiple seeds model
resultsOutSQMs2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs2$OOB_RMSE
[1] 5.532929
resultsOutSQMs2$OOB_MAE
[1] 4.391987
resultsOutSQMs2$Rsquared
[1] 0.6285529
# store results
resdHiAnnoyFit['Abs SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resdHiAnnoyFit['Abs SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_MAE
resdHiAnnoyFit['Abs SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resdHiAnnoyPermImp$AbsSQMs2 <- resultsOutSQMs2$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 50))
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyAbsSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyAbsSQMsNoTonLdConPermimp.svg")
ggsave(filename="PtAdHiAnnoyAbsSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyAbsSQMsNoTonLdConPermimp.pdf")
}
Next, the difference metrics are analysed,
Here, the ‘ambient only’ stimuli are removed, as the analysis cannot handle missing values for dB metrics.
stimDataANum <- stimDataANum[complete.cases(stimDataANum),]
stimDataANum$UASLAeq <- as.numeric(stimDataANum$UASLAeq)
stimDataANum$SNRlevel <- as.numeric(stimDataANum$SNRlevel)
iVars <- c(names(stimDataANum)[which(colnames(stimDataANum)=="LAeqLAF90diff"):
which(colnames(stimDataANum)=="dImpulsSHM05ExMaxLR")], "SNRlevel")
dVar <- "dHighAnnoyPc"
seeds <- c(568392, 498, 4089, 78132, 741809)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 4001
mtry <- as.integer(length(iVars)/2.25)
Train preliminary model
nperm <- 30
resultsOutDiffs <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutDiffs$OOB_RMSE
[1] 5.719733
resultsOutDiffs$OOB_MAE
[1] 4.428596
resultsOutDiffs$Rsquared
[1] 0.6091806
Train multiple seeds model
resultsOutDiffs <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutDiffs$OOB_RMSE
[1] 5.715133
resultsOutDiffs$OOB_MAE
[1] 4.425906
resultsOutDiffs$Rsquared
[1] 0.609705
# store results
resdHiAnnoyFit['Diff vars', 'RMSE'] <- resultsOutDiffs$OOB_RMSE
resdHiAnnoyFit['Diff vars', 'MAE'] <- resultsOutDiffs$OOB_MAE
resdHiAnnoyFit['Diff vars', 'Rsquared'] <- resultsOutDiffs$Rsquared
resdHiAnnoyPermImp$DiffVars <- resultsOutDiffs$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutDiffs.conimp <- arrange(resultsOutDiffs$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimp), levels=rownames(resultsOutDiffs.conimp)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyDiffVarsConPermimp.svg", width=8, height=10, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyDiffVarsConPermimp.svg")
ggsave(filename="PtAdHiAnnoyDiffVarsConPermimp.pdf", width=8, height=10, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyDiffVarsConPermimp.pdf")
}
# Plot only positive values
resultsOutDiffs.conimpPtv <- resultsOutDiffs.conimp |>
rownames_to_column('Metric') |>
filter_if(is.numeric, all_vars(. > 0)) |>
column_to_rownames('Metric')
pBar <- ggplot(resultsOutDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimpPtv), levels=rownames(resultsOutDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyDiffVarsConPermimpPtv.svg", width=8, height=8, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyDiffVarsConPermimpPtv.svg")
ggsave(filename="PtAdHiAnnoyDiffVarsConPermimpPtv.pdf", width=8, height=8, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyDiffVarsConPermimpPtv.pdf")
}
Selected metric
diffVar <- "EPNLLAeqdiff"
The ambient-only stimuli are NOT replaced here, as the difference analysis includes dB metrics.
iVars <- c(diffVar, "dSharpAurISO3PowAvgMaxLR", "dSharpAurISO305ExMaxLR", "dSharpAurSHMPowAvgMaxLR", "dSharpAurSHM05ExMaxLR")
dVar <- "dHighAnnoyPc"
seeds <- c(84194, 905, 64815, 928054, 625091, 582031)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutSharp <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSharp$OOB_RMSE
[1] 5.817107
resultsOutSharp$OOB_MAE
[1] 4.431913
resultsOutSharp$Rsquared
[1] 0.5844272
Train multiple seeds model
resultsOutSharp <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSharp$OOB_RMSE
[1] 5.794601
resultsOutSharp$OOB_MAE
[1] 4.410331
resultsOutSharp$Rsquared
[1] 0.5874878
# store results
resdHiAnnoyFit['Diff sharp', 'RMSE'] <- resultsOutSharp$OOB_RMSE
resdHiAnnoyFit['Diff sharp', 'MAE'] <- resultsOutSharp$OOB_MAE
resdHiAnnoyFit['Diff sharp', 'Rsquared'] <- resultsOutSharp$Rsquared
resdHiAnnoyPermImp$DiffSharp <- resultsOutSharp$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSharp.conimp <- arrange(resultsOutSharp$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSharp.conimp) + geom_col(aes(x=factor(rownames(resultsOutSharp.conimp), levels=rownames(resultsOutSharp.conimp)), y=CondPermImp), fill=mycolours[2], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dSharpness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoydSharpConPermimp.svg", width=8, height=2.6, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoydSharpConPermimp.svg")
ggsave(filename="PtAdHiAnnoydSharpConPermimp.pdf", width=8, height=2.6, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoydSharpConPermimp.pdf")
}
Selected metric
dSharpVar <- "dSharpAurSHMPowAvgMaxLR"
iVars <- c(diffVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR", "dTonLdECMAPowAvgBin", "dTonLdECMA05ExBin")
dVar <- "dHighAnnoyPc"
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
seeds <- c(561684, 104798, 1536, 48, 48561)
Selected hyperparameters
ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
# Tonality with tonal loudness
nperm <- 5
resultsOutTonal1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal1$OOB_RMSE
[1] 6.71205
resultsOutTonal1$OOB_MAE
[1] 5.189994
resultsOutTonal1$Rsquared
[1] 0.4460736
Train multiple seeds model
# Tonality with tonal loudness
resultsOutTonal1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal1$OOB_RMSE
[1] 6.722309
resultsOutTonal1$OOB_MAE
[1] 5.199688
resultsOutTonal1$Rsquared
[1] 0.4443759
# store results
resdHiAnnoyFit['Diff tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resdHiAnnoyFit['Diff tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resdHiAnnoyFit['Diff tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resdHiAnnoyPermImp$DiffTonal1 <- resultsOutTonal1$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dTonality inc. dtonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 70))
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoydTonalLdConPermimp.svg", width=8, height=2.6, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoydTonalLdConPermimp.svg")
ggsave(filename="PtAdHiAnnoydTonalLdConPermimp.pdf", width=8, height=2.6, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoydTonalLdConPermimp.pdf")
}
Selected metric
dTonLdVar <- "dTonLdECMAPowAvgBin"
iVars <- c(diffVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR")
dVar <- "dHighAnnoyPc"
seeds <- c(410865, 2954, 70812, 203, 7984)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
NA
NA
Selected hyperparameters
ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
# Tonality
nperm <- 5
resultsOutTonal2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal2$OOB_RMSE
[1] 6.581687
resultsOutTonal2$OOB_MAE
[1] 5.067818
resultsOutTonal2$Rsquared
[1] 0.4671435
Train multiple seeds model
# Tonality
resultsOutTonal2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal2$OOB_RMSE
[1] 6.580352
resultsOutTonal2$OOB_MAE
[1] 5.068529
resultsOutTonal2$Rsquared
[1] 0.4673965
# store results
resdHiAnnoyFit['Diff tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resdHiAnnoyFit['Diff tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resdHiAnnoyFit['Diff tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resdHiAnnoyPermImp$DiffTonal2 <- resultsOutTonal2$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dTonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 70))
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoydTonalConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoydTonalConPermimp.svg")
ggsave(filename="PtAdHiAnnoydTonalConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoydTonalConPermimp.pdf")
}
Selected metric
dTonalVar <- "dTonalSHMIntAvgMaxLR"
# Fluctuation strength
iVars <- c(diffVar, "dFluctSHM10ExBin", "dFluctSHM05ExBin", "dFluctOV10ExMaxLR", "dFluctOV05ExMaxLR")
dVar <- "dHighAnnoyPc"
seeds <- c(418657, 84, 1630, 18659, 3687)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 251
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutFluct <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutFluct$OOB_RMSE
[1] 6.714071
resultsOutFluct$OOB_MAE
[1] 5.1887
resultsOutFluct$Rsquared
[1] 0.4455587
Train multiple seeds model
resultsOutFluct <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutFluct$OOB_RMSE
[1] 6.702409
resultsOutFluct$OOB_MAE
[1] 5.191661
resultsOutFluct$Rsquared
[1] 0.4472641
# store results
resdHiAnnoyFit['Diff fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resdHiAnnoyFit['Diff fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resdHiAnnoyFit['Diff fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resdHiAnnoyPermImp$DiffFluct <- resultsOutFluct$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dFluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoydFluctConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyFluctConPermimp.svg")
ggsave(filename="PtAdHiAnnoydFluctConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyFluctConPermimp.pdf")
}
Selected metric
dFluctVar <- "dFluctSHM10ExBin"
# Roughness
iVars <- c(diffVar, "dRoughECMA10ExBin", "dRoughECMA05ExBin", "dRoughFZ10ExMaxLR", "dRoughFZ05ExMaxLR")
dVar <- "dHighAnnoyPc"
seeds <- c(69851, 85109, 410986, 1563, 896)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutRough <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
Warning: Permutation importance rank order within 10% of max differs between seeds: increase number of trees ('ntree') or number of permutations ('nperm'), or subsample of features ('mtry')
# print model prediction results
resultsOutRough$OOB_RMSE
NULL
resultsOutRough$OOB_MAE
NULL
resultsOutRough$Rsquared
NULL
Train multiple seeds model
resultsOutRough <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutRough$OOB_RMSE
[1] 6.950843
resultsOutRough$OOB_MAE
[1] 5.259601
resultsOutRough$Rsquared
[1] 0.4054285
# store results
resdHiAnnoyFit['Diff rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resdHiAnnoyFit['Diff rough', 'MAE'] <- resultsOutRough$OOB_MAE
resdHiAnnoyFit['Diff rough', 'Rsquared'] <- resultsOutRough$Rsquared
resdHiAnnoyPermImp$DiffRough <- resultsOutRough$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dRoughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoydRoughConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoydRoughConPermimp.svg")
ggsave(filename="PtAdHiAnnoydRoughConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoydRoughConPermimp.pdf")
}
Selected metric
dRoughVar <- "dRoughFZ05ExMaxLR"
# Impulsiveness
iVars <- c(diffVar, "dImpulsSHMAvgMaxLR", "dImpulsSHM05ExMaxLR", "dImpulsSHMPowAvgMaxLR")
dVar <- "dHighAnnoyPc"
seeds <- c(418659, 7805, 38475, 65834, 1653)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutImpuls <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutImpuls$OOB_RMSE
[1] 6.674278
resultsOutImpuls$OOB_MAE
[1] 5.175702
resultsOutImpuls$Rsquared
[1] 0.4517669
Train multiple seeds model
resultsOutImpuls <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutImpuls$OOB_RMSE
[1] 6.673098
resultsOutImpuls$OOB_MAE
[1] 5.185144
resultsOutImpuls$Rsquared
[1] 0.4520846
# store results
resdHiAnnoyFit['Diff impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resdHiAnnoyFit['Diff impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resdHiAnnoyFit['Diff impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resdHiAnnoyPermImp$DiffImpuls <- resultsOutImpuls$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + ggtitle("dImpulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoydImpulsConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoydImpulsConPermimp.svg")
ggsave(filename="PtAdHiAnnoydImpulsConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoydImpulsConPermimp.pdf")
}
Selected metric
dImpulsVar <- "dImpulsSHMAvgMaxLR"
Now the highest importance dSQMs are ranked against each other, controlling for loudness difference.
iVars <- c(diffVar, dSharpVar, dTonLdVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "dHighAnnoyPc"
seeds <- c(98465, 54163, 6541, 36485, 849675)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 2501
mtry <- as.integer(length(iVars)/1.75)
Train preliminary model
nperm <- 5
resultsOutSQMs1 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
Warning: Permutation importance rank order within 10% of max differs between seeds: increase number of trees ('ntree') or number of permutations ('nperm'), or subsample of features ('mtry')
# print model prediction results
resultsOutSQMs1$OOB_RMSE
NULL
resultsOutSQMs1$OOB_MAE
NULL
resultsOutSQMs1$Rsquared
NULL
Train multiple seeds model
resultsOutSQMs1 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs1$OOB_RMSE
[1] 5.939049
resultsOutSQMs1$OOB_MAE
[1] 4.654886
resultsOutSQMs1$Rsquared
[1] 0.5768752
# store results
resdHiAnnoyFit['Diff SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resdHiAnnoyFit['Diff SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resdHiAnnoyFit['Diff SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resdHiAnnoyPermImp$DiffSQMs1 <- resultsOutSQMs1$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 50))
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyDiffSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyDiffSQMsTonLdConPermimp.svg")
ggsave(filename="PtAdHiAnnoyDiffSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyDiffSQMsTonLdConPermimp.pdf")
}
iVars <- c(diffVar, dSharpVar, dTonalVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "dHighAnnoyPc"
seeds <- c(49865, 7852, 845961, 410583, 36748)
p <- mtryTune(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 4001
mtry <- as.integer(length(iVars)/1.75)
Train preliminary model
nperm <- 5
resultsOutSQMs2 <- crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs2$OOB_RMSE
[1] 5.955857
resultsOutSQMs2$OOB_MAE
[1] 4.656295
resultsOutSQMs2$Rsquared
[1] 0.5742314
Train multiple seeds model
resultsOutSQMs2 <- multi_crfReg(dataIn=stimDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs2$OOB_RMSE
[1] 5.966773
resultsOutSQMs2$OOB_MAE
[1] 4.669058
resultsOutSQMs2$Rsquared
[1] 0.5722321
# store results
resdHiAnnoyFit['Diff SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resdHiAnnoyFit['Diff SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_MAE
resdHiAnnoyFit['Diff SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resdHiAnnoyPermImp$DiffSQMs2 <- resultsOutSQMs2$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% highly annoyed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 50))
pBar
if (saveplots){
ggsave(filename="PtAdHiAnnoyDiffSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
unlink("PtAdHiAnnoyDiffSQMsNoTonLdConPermimp.svg")
ggsave(filename="PtAdHiAnnoyDiffSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
unlink("PtAdHiAnnoyDiffSQMsNoTonLdConPermimp.pdf")
}
if (savedata){
utils::write.csv(resdHiAnnoyFit, paste(outDataPath, "\\ptACRFdHiAnnoyOOBFit.csv", sep=""))
ii <- 0
temp = list()
for (res in resdHiAnnoyPermImp){
ii <- ii + 1
temp[[ii]] <- as.data.frame(resdHiAnnoyPermImp[ii])
names(temp[[ii]]) <- names(resdHiAnnoyPermImp[ii])
}
openxlsx::write.xlsx(temp, paste(outDataPath, "\\ptACRFdHiAnnoyConPermimp.xlsx",
sep=""),
rowNames=TRUE)
}
Initialise results output variables
resNoticeFit <- data.frame(RMSE = numeric(),
MAE = numeric(),
Rsquared = numeric())
resNoticePermImp <- list()
Here, the ‘ambient only’ stimuli are removed, as the analysis cannot handle missing values for dB metrics.
stimNoticeDataANum <- stimNoticeDataANum[complete.cases(stimNoticeDataANum),]
stimNoticeDataANum$UASLAeq <- as.numeric(stimNoticeDataANum$UASLAeq)
stimNoticeDataANum$SNRlevel <- as.numeric(stimNoticeDataANum$SNRlevel)
iVars <- names(stimNoticeDataANum)[which(names(stimNoticeDataANum) == 'UASLAeq'):which(names(stimNoticeDataANum) == 'UASImpulsSHM05ExMaxLR')]
iVars <- iVars[! iVars %in% 'SNRlevel']
iVars <- c(iVars,
names(stimNoticeDataANum)[which(colnames(stimNoticeDataANum)=="LAeqLAF90diff"):
which(colnames(stimNoticeDataANum)=="dImpulsSHM05ExMaxLR")], "SNRlevel")
dVar <- "NoticedPcFilt"
seeds <- c(2, 312, 1897, 465978, 821659)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 1501
mtry <- as.integer(length(iVars)/2.75)
Train preliminary model
nperm <- 10
resultsOutAbsDiffs <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutAbsDiffs$OOB_RMSE
[1] 10.23042
resultsOutAbsDiffs$OOB_MAE
[1] 7.618346
resultsOutAbsDiffs$Rsquared
[1] 0.9097534
Train multiple seeds model
resultsOutAbsDiffs <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutAbsDiffs$OOB_RMSE
[1] 10.19738
resultsOutAbsDiffs$OOB_MAE
[1] 7.652439
resultsOutAbsDiffs$Rsquared
[1] 0.9103752
# store results
resNoticeFit['All vars', 'RMSE'] <- resultsOutAbsDiffs$OOB_RMSE
resNoticeFit['All vars', 'MAE'] <- resultsOutAbsDiffs$OOB_MAE
resNoticeFit['All vars', 'Rsquared'] <- resultsOutAbsDiffs$Rsquared
resNoticePermImp$AllVars <- resultsOutAbsDiffs$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutAbsDiffs.conimp <- arrange(resultsOutAbsDiffs$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutAbsDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimp), levels=rownames(resultsOutAbsDiffs.conimp)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtANoticeAllVarsConPermimp.svg", width=8, height=20, path=file.path(outFigPath, "svg"))
unlink("PtANoticeAllVarsConPermimp.svg")
ggsave(filename="PtANoticeAllVarsConPermimp.pdf", width=8, height=20, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeAllVarsConPermimp.pdf")
}
# Plot only positive values
resultsOutAbsDiffs.conimpPtv <- resultsOutAbsDiffs.conimp |>
rownames_to_column('Metric') |>
filter_if(is.numeric, all_vars(. > 0)) |>
column_to_rownames('Metric')
pBar <- ggplot(resultsOutAbsDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutAbsDiffs.conimpPtv), levels=rownames(resultsOutAbsDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[9], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtANoticeAllVarsConPermimpPtv.svg", width=8, height=12, path=file.path(outFigPath, "svg"))
unlink("PtANoticeAllVarsConPermimp.svg")
ggsave(filename="PtANoticeAllVarsConPermimpPtv.pdf", width=8, height=12, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeAllVarsConPermimp.pdf")
}
iVars <- names(stimNoticeDataANum)[which(names(stimNoticeDataANum) == 'UASLAeq'):which(names(stimNoticeDataANum) == 'UASImpulsSHM05ExMaxLR')]
iVars <- iVars[! iVars %in% c('SNRlevel')]
dVar <- "NoticedPcFilt"
seeds <- c(578312, 544, 84894, 54654, 153157)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
if (saveplots){
ggsave(filename="PtANoticeAbsVarsHyperTune.svg", width=12, height=4, path=file.path(outFigPath, "svg"))
unlink("PtANoticeAbsVarsHyperTune.svg")
ggsave(filename="PtANoticeAbsVarsHyperTune.pdf", width=12, height=4, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeAbsVarsHyperTune.pdf")
}
Selected hyperparameters
ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 20
resultsOutAbs <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutAbs$OOB_RMSE
[1] 16.59018
resultsOutAbs$OOB_MAE
[1] 12.64339
resultsOutAbs$Rsquared
[1] 0.7829788
Train multiple seeds model
resultsOutAbs <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutAbs$OOB_RMSE
[1] 16.50718
resultsOutAbs$OOB_MAE
[1] 12.55933
resultsOutAbs$Rsquared
[1] 0.7846705
# store results
resNoticeFit['Abs vars', 'RMSE'] <- resultsOutAbs$OOB_RMSE
resNoticeFit['Abs vars', 'MAE'] <- resultsOutAbs$OOB_MAE
resNoticeFit['Abs vars', 'Rsquared'] <- resultsOutAbs$Rsquared
resNoticePermImp$AbsVars <- resultsOutAbs$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutAbs.conimp <- arrange(resultsOutAbs$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutAbs.conimp) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimp), levels=rownames(resultsOutAbs.conimp)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) +
coord_flip()
pBar
if (saveplots){
ggsave(filename="PtANoticeAbsVarsConPermimp.svg", width=8, height=11.5, path=file.path(outFigPath, "svg"))
unlink("PtANoticeAbsVarsConPermimp.svg")
ggsave(filename="PtANoticeAbsVarsConPermimp.pdf", width=8, height=11.5, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeAbsVarsConPermimp.pdf")
}
# Plot only positive values
resultsOutAbs.conimpPtv <- resultsOutAbs.conimp |>
rownames_to_column('Metric') |>
filter_if(is.numeric, all_vars(. > 0)) |>
column_to_rownames('Metric')
pBar <- ggplot(resultsOutAbs.conimpPtv,) + geom_col(aes(x=factor(rownames(resultsOutAbs.conimpPtv), levels=rownames(resultsOutAbs.conimpPtv)), y=CondPermImp), fill=mycolours[1], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtANoticeAbsVarsConPermimpPtv.svg", width=8, height=7, path=file.path(outFigPath, "svg"))
unlink("PtANoticeAbsVarsConPermimpPtv.svg")
ggsave(filename="PtANoticeAbsVarsConPermimpPtv.pdf", width=8, height=7, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeAbsVarsConPermimpPtv.pdf")
}
Selected metric
absVar <- "UASPNLmaxMaxLR"
Here, the ‘ambient only’ stimuli are not replaced, as the analysis cannot handle missing values for dB metrics (PNL).
#
# stimNoticeDataANum <- rbind(stimNoticeDataANum, stimNoticeDataA[stimNoticeDataA['Stimulus']
# == "A1",
# colnames(stimNoticeDataANum)],
# stimNoticeDataA[stimNoticeDataA['Stimulus']
# == "A2",
# colnames(stimNoticeDataANum)])
# stimNoticeDataANum <- arrange(stimNoticeDataANum, Stimulus)
iVars <- c(absVar, "AmbientLAeq", "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonLdECMAPowAvgBin", "UASTonLdECMA05ExBin", "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR")
dVar <- "NoticedPcFilt"
seeds <- c(540, 104798, 456464, 87331, 94564)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
# Tonality with tonal loudness
nperm <- 5
resultsOutTonal1 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal1$OOB_RMSE
[1] 18.47209
resultsOutTonal1$OOB_MAE
[1] 14.20136
resultsOutTonal1$Rsquared
[1] 0.7179216
Train multiple seeds model
# Tonality with tonal loudness
resultsOutTonal1 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal1$OOB_RMSE
[1] 18.19492
resultsOutTonal1$OOB_MAE
[1] 14.0021
resultsOutTonal1$Rsquared
[1] 0.7280624
# store results
resNoticeFit['Abs tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resNoticeFit['Abs tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resNoticeFit['Abs tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resNoticePermImp$AbsTonal1 <- resultsOutTonal1$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("Tonality inc. tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 1000))
pBar
if (saveplots){
ggsave(filename="PtANoticeTonalLdConPermimp.svg", width=8, height=3.8, path=file.path(outFigPath, "svg"))
unlink("PtANoticeTonalLdConPermimp.svg")
ggsave(filename="PtANoticeTonalLdConPermimp.pdf", width=8, height=3.8, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeTonalLdConPermimp.pdf")
}
Selected metric
tonLdVar <- "UASTonLdECMA05ExBin"
iVars <- c(absVar, "AmbientLAeq", "UASTonalECMAAvgMaxLR", "UASTonalSHMInt05ExMaxLR", "UASTonalSHMIntAvgMaxLR", "UASTonalECMA05ExMaxLR", "UASTonalAurAvgMaxLR", "UASTonalAur05ExMaxLR", "UASTonalAur10ExMaxLR")
dVar <- "NoticedPcFilt"
seeds <- c(156089, 5860, 10528, 89541, 4685146)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
# Tonality
nperm <- 5
resultsOutTonal2 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal2$OOB_RMSE
[1] 17.0126
resultsOutTonal2$OOB_MAE
[1] 13.20926
resultsOutTonal2$Rsquared
[1] 0.7729116
Train multiple seeds model
# Tonality
resultsOutTonal2 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal2$OOB_RMSE
[1] 17.06205
resultsOutTonal2$OOB_MAE
[1] 13.236
resultsOutTonal2$Rsquared
[1] 0.7718614
# store results
resNoticeFit['Abs tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resNoticeFit['Abs tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resNoticeFit['Abs tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resNoticePermImp$AbsTonal2 <- resultsOutTonal2$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("Tonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 1000))
pBar
if (saveplots){
ggsave(filename="PtANoticeTonalConPermimp.svg", width=8, height=3.2, path=file.path(outFigPath, "svg"))
unlink("PtANoticeTonalConPermimp.svg")
ggsave(filename="PtANoticeTonalConPermimp.pdf", width=8, height=3.2, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeTonalConPermimp.pdf")
}
Selected metric
tonalVar <- "UASTonalSHMInt05ExMaxLR"
# Fluctuation strength
iVars <- c(absVar, "AmbientLAeq", "UASFluctSHM10ExBin", "UASFluctSHM05ExBin", "UASFluctFZ10ExMaxLR", "UASFluctFZ05ExMaxLR", "UASFluctOV10ExMaxLR", "UASFluctOV05ExMaxLR")
dVar <- "NoticedPcFilt"
seeds <- c(25107, 546098, 195, 5937, 102658)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 251
mtry <- as.integer(length(iVars)/1.5)
Train preliminary model
nperm <- 5
resultsOutFluct <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutFluct$OOB_RMSE
[1] 16.91354
resultsOutFluct$OOB_MAE
[1] 12.93498
resultsOutFluct$Rsquared
[1] 0.7866218
Train multiple seeds model
resultsOutFluct <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutFluct$OOB_RMSE
[1] 16.96583
resultsOutFluct$OOB_MAE
[1] 12.82457
resultsOutFluct$Rsquared
[1] 0.7842884
# store results
resNoticeFit['Abs fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resNoticeFit['Abs fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resNoticeFit['Abs fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resNoticePermImp$AbsFluct <- resultsOutFluct$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("Fluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtANoticeFluctConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
unlink("PtANoticeFluctConPermimp.svg")
ggsave(filename="PtANoticeFluctConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeFluctConPermimp.pdf")
}
Selected metric
fluctVar <- "UASFluctOV10ExMaxLR"
# Roughness
iVars <- c(absVar, "AmbientLAeq", "UASRoughECMA10ExBin", "UASRoughECMA05ExBin", "UASRoughFZ10ExMaxLR", "UASRoughFZ05ExMaxLR", "UASRoughDW10ExMaxLR", "UASRoughDW05ExMaxLR")
dVar <- "NoticedPcFilt"
seeds <- c(4701, 52187, 16589, 65217, 16893)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 1001
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutRough <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutRough$OOB_RMSE
[1] 17.26411
resultsOutRough$OOB_MAE
[1] 13.17981
resultsOutRough$Rsquared
[1] 0.7714465
Train multiple seeds model
resultsOutRough <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutRough$OOB_RMSE
[1] 16.75907
resultsOutRough$OOB_MAE
[1] 12.82419
resultsOutRough$Rsquared
[1] 0.7847338
# store results
resNoticeFit['Abs rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resNoticeFit['Abs rough', 'MAE'] <- resultsOutRough$OOB_MAE
resNoticeFit['Abs rough', 'Rsquared'] <- resultsOutRough$Rsquared
resNoticePermImp$AbsRough <- resultsOutRough$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("Roughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtANoticeRoughConPermimp.svg", width=8, height=2.9, path=file.path(outFigPath, "svg"))
unlink("PtANoticeRoughConPermimp.svg")
ggsave(filename="PtANoticeRoughConPermimp.pdf", width=8, height=2.9, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeRoughConPermimp.pdf")
}
Selected metric
roughVar <- "UASRoughFZ05ExMaxLR"
# Impulsiveness
iVars <- c(absVar, "AmbientLAeq", "UASImpulsSHMAvgMaxLR", "UASImpulsSHM05ExMaxLR", "UASImpulsSHMPowAvgMaxLR")
dVar <- "NoticedPcFilt"
seeds <- c(8495, 59867, 5416, 9843, 86)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 5501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutImpuls <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutImpuls$OOB_RMSE
[1] 16.24906
resultsOutImpuls$OOB_MAE
[1] 12.47619
resultsOutImpuls$Rsquared
[1] 0.7966846
Train multiple seeds model
resultsOutImpuls <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutImpuls$OOB_RMSE
[1] 16.26418
resultsOutImpuls$OOB_MAE
[1] 12.46156
resultsOutImpuls$Rsquared
[1] 0.7959352
# store results
resNoticeFit['Abs impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resNoticeFit['Abs impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resNoticeFit['Abs impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resNoticePermImp$AbsImpuls <- resultsOutImpuls$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("Impulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtANoticeImpulsConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtANoticeImpulsConPermimp.svg")
ggsave(filename="PtANoticeImpulsConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeImpulsConPermimp.pdf")
}
Selected metric
impulsVar <- "UASImpulsSHMAvgMaxLR"
Now the highest importance SQMs are ranked against each other, controlling for UAS loudness and ambient LAeq.
iVars <- c(absVar, "AmbientLAeq", sharpVar, tonLdVar, fluctVar, roughVar, impulsVar)
dVar <- "NoticedPcFilt"
seeds <- c(70498, 4, 14986, 453, 864)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 251
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutSQMs1 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs1$OOB_RMSE
[1] 18.88386
resultsOutSQMs1$OOB_MAE
[1] 14.50936
resultsOutSQMs1$Rsquared
[1] 0.7014957
Train multiple seeds model
resultsOutSQMs1 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs1$OOB_RMSE
[1] 18.76572
resultsOutSQMs1$OOB_MAE
[1] 14.41058
resultsOutSQMs1$Rsquared
[1] 0.7076868
# store results
resNoticeFit['Abs SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resNoticeFit['Abs SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resNoticeFit['Abs SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resNoticePermImp$AbsSQMs1 <- resultsOutSQMs1$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 800))
pBar
if (saveplots){
ggsave(filename="PtANoticeAbsSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
unlink("PtANoticeAbsSQMsTonLdConPermimp.svg")
ggsave(filename="PtANoticeAbsSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeAbsSQMsTonLdConPermimp.pdf")
}
iVars <- c(absVar, "AmbientLAeq", sharpVar, tonalVar, fluctVar, roughVar, impulsVar)
dVar <- "NoticedPcFilt"
seeds <- c(546, 57203, 270835, 60592, 8094)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutSQMs2 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs2$OOB_RMSE
[1] 18.88129
resultsOutSQMs2$OOB_MAE
[1] 14.46072
resultsOutSQMs2$Rsquared
[1] 0.7037345
Train multiple seeds model
resultsOutSQMs2 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs2$OOB_RMSE
[1] 18.78148
resultsOutSQMs2$OOB_MAE
[1] 14.40389
resultsOutSQMs2$Rsquared
[1] 0.7078282
# store results
resNoticeFit['Abs SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resNoticeFit['Abs SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_MAE
resNoticeFit['Abs SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resNoticePermImp$AbsSQMs2 <- resultsOutSQMs2$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(-1, 800))
pBar
if (saveplots){
ggsave(filename="PtANoticeAbsSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
unlink("PtANoticeAbsSQMsNoTonLdConPermimp.svg")
ggsave(filename="PtANoticeAbsSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeAbsSQMsNoTonLdConPermimp.pdf")
}
Next, the difference metrics are analysed,
Here, the ‘ambient only’ stimuli are removed, as the analysis cannot handle missing values for dB metrics.
stimNoticeDataANum <- stimNoticeDataANum[complete.cases(stimNoticeDataANum),]
stimNoticeDataANum$UASLAeq <- as.numeric(stimNoticeDataANum$UASLAeq)
stimNoticeDataANum$SNRlevel <- as.numeric(stimNoticeDataANum$SNRlevel)
iVars <- c(names(stimNoticeDataANum)[which(colnames(stimNoticeDataANum)=="LAeqLAF90diff"):
which(colnames(stimNoticeDataANum)=="dImpulsSHM05ExMaxLR")], "SNRlevel")
dVar <- "NoticedPcFilt"
seeds <- c(568392, 498, 4089, 78132, 741809)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 1501
mtry <- as.integer(length(iVars)/3.5)
Train preliminary model
nperm <- 30
resultsOutDiffs <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutDiffs$OOB_RMSE
[1] 10.09164
resultsOutDiffs$OOB_MAE
[1] 7.746637
resultsOutDiffs$Rsquared
[1] 0.9128949
Train multiple seeds model
resultsOutDiffs <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutDiffs$OOB_RMSE
[1] 10.15461
resultsOutDiffs$OOB_MAE
[1] 7.781117
resultsOutDiffs$Rsquared
[1] 0.9118708
# store results
resNoticeFit['Diff vars', 'RMSE'] <- resultsOutDiffs$OOB_RMSE
resNoticeFit['Diff vars', 'MAE'] <- resultsOutDiffs$OOB_MAE
resNoticeFit['Diff vars', 'Rsquared'] <- resultsOutDiffs$Rsquared
resNoticePermImp$DiffVars <- resultsOutDiffs$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutDiffs.conimp <- arrange(resultsOutDiffs$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutDiffs.conimp) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimp), levels=rownames(resultsOutDiffs.conimp)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtANoticeDiffVarsConPermimp.svg", width=8, height=10, path=file.path(outFigPath, "svg"))
unlink("PtANoticeDiffVarsConPermimp.svg")
ggsave(filename="PtANoticeDiffVarsConPermimp.pdf", width=8, height=10, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeDiffVarsConPermimp.pdf")
}
# Plot only positive values
resultsOutDiffs.conimpPtv <- resultsOutDiffs.conimp |>
rownames_to_column('Metric') |>
filter_if(is.numeric, all_vars(. > 0)) |>
column_to_rownames('Metric')
pBar <- ggplot(resultsOutDiffs.conimpPtv) + geom_col(aes(x=factor(rownames(resultsOutDiffs.conimpPtv), levels=rownames(resultsOutDiffs.conimpPtv)), y=CondPermImp), fill=mycolours[8], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtANoticeDiffVarsConPermimpPtv.svg", width=8, height=8, path=file.path(outFigPath, "svg"))
unlink("PtANoticeDiffVarsConPermimpPtv.svg")
ggsave(filename="PtANoticeDiffVarsConPermimpPtv.pdf", width=8, height=8, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeDiffVarsConPermimpPtv.pdf")
}
Selected metric
diffVar <- "Detect0p5dBMaxMaxLR"
The ambient-only stimuli are NOT replaced here, as the difference analysis includes dB metrics.
iVars <- c(diffVar, "dSharpAurISO3PowAvgMaxLR", "dSharpAurISO305ExMaxLR", "dSharpAurSHMPowAvgMaxLR", "dSharpAurSHM05ExMaxLR")
dVar <- "NoticedPcFilt"
seeds <- c(84194, 905, 928054, 625091, 582031)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutSharp <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSharp$OOB_RMSE
[1] 10.40193
resultsOutSharp$OOB_MAE
[1] 8.108278
resultsOutSharp$Rsquared
[1] 0.9042801
Train multiple seeds model
resultsOutSharp <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSharp$OOB_RMSE
[1] 10.43319
resultsOutSharp$OOB_MAE
[1] 8.137405
resultsOutSharp$Rsquared
[1] 0.9035948
# store results
resNoticeFit['Diff sharp', 'RMSE'] <- resultsOutSharp$OOB_RMSE
resNoticeFit['Diff sharp', 'MAE'] <- resultsOutSharp$OOB_MAE
resNoticeFit['Diff sharp', 'Rsquared'] <- resultsOutSharp$Rsquared
resNoticePermImp$DiffSharp <- resultsOutSharp$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSharp.conimp <- arrange(resultsOutSharp$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSharp.conimp) + geom_col(aes(x=factor(rownames(resultsOutSharp.conimp), levels=rownames(resultsOutSharp.conimp)), y=CondPermImp), fill=mycolours[2], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("dSharpness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtANoticedSharpConPermimp.svg", width=8, height=2.6, path=file.path(outFigPath, "svg"))
unlink("PtANoticedSharpConPermimp.svg")
ggsave(filename="PtANoticedSharpConPermimp.pdf", width=8, height=2.6, path=file.path(outFigPath, "pdf"))
unlink("PtANoticedSharpConPermimp.pdf")
}
Selected metric
sharpVar <- "dSharpAurISO305ExMaxLR"
iVars <- c(diffVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR", "dTonLdECMAPowAvgBin", "dTonLdECMA05ExBin")
dVar <- "NoticedPcFilt"
seeds <- c(561684, 104798, 1536, 48, 48561)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <-1501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
# Tonality with tonal loudness
nperm <- 5
resultsOutTonal1 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal1$OOB_RMSE
[1] 10.46408
resultsOutTonal1$OOB_MAE
[1] 8.013044
resultsOutTonal1$Rsquared
[1] 0.9021881
Train multiple seeds model
# Tonality with tonal loudness
resultsOutTonal1 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal1$OOB_RMSE
[1] 10.45256
resultsOutTonal1$OOB_MAE
[1] 8.001321
resultsOutTonal1$Rsquared
[1] 0.9026108
# store results
resNoticeFit['Diff tonal inc loud', 'RMSE'] <- resultsOutTonal1$OOB_RMSE
resNoticeFit['Diff tonal inc loud', 'MAE'] <- resultsOutTonal1$OOB_MAE
resNoticeFit['Diff tonal inc loud', 'Rsquared'] <- resultsOutTonal1$Rsquared
resNoticePermImp$DiffTonal1 <- resultsOutTonal1$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutTonal1.conimp <- arrange(resultsOutTonal1$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutTonal1.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal1.conimp), levels=rownames(resultsOutTonal1.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("dTonality inc. dtonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1200))
pBar
if (saveplots){
ggsave(filename="PtANoticedTonalLdConPermimp.svg", width=8, height=2.6, path=file.path(outFigPath, "svg"))
unlink("PtANoticedTonalLdConPermimp.svg")
ggsave(filename="PtANoticedTonalLdConPermimp.pdf", width=8, height=2.6, path=file.path(outFigPath, "pdf"))
unlink("PtANoticedTonalLdConPermimp.pdf")
}
Selected metric
tonLdVar <- "dTonLdECMAPowAvgBin"
iVars <- c(diffVar, "dTonalECMAAvgMaxLR", "dTonalSHMInt05ExMaxLR", "dTonalSHMIntAvgMaxLR", "dTonalECMA05ExMaxLR")
dVar <- "NoticedPcFilt"
seeds <- c(410865, 2954, 70812, 203, 7984)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
# Tonality
nperm <- 5
resultsOutTonal2 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal2$OOB_RMSE
[1] 12.53815
resultsOutTonal2$OOB_MAE
[1] 9.798852
resultsOutTonal2$Rsquared
[1] 0.8608015
Train multiple seeds model
# Tonality
resultsOutTonal2 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutTonal2$OOB_RMSE
[1] 12.57501
resultsOutTonal2$OOB_MAE
[1] 9.807799
resultsOutTonal2$Rsquared
[1] 0.8595653
# store results
resNoticeFit['Diff tonal no loud', 'RMSE'] <- resultsOutTonal2$OOB_RMSE
resNoticeFit['Diff tonal no loud', 'MAE'] <- resultsOutTonal2$OOB_MAE
resNoticeFit['Diff tonal no loud', 'Rsquared'] <- resultsOutTonal2$Rsquared
resNoticePermImp$DiffTonal2 <- resultsOutTonal2$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutTonal2.conimp <- arrange(resultsOutTonal2$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutTonal2.conimp) + geom_col(aes(x=factor(rownames(resultsOutTonal2.conimp), levels=rownames(resultsOutTonal2.conimp)), y=CondPermImp), fill=mycolours[3], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("dTonality w/o tonal loudness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 1200))
pBar
if (saveplots){
ggsave(filename="PtANoticedTonalConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtANoticedTonalConPermimp.svg")
ggsave(filename="PtANoticedTonalConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtANoticedTonalConPermimp.pdf")
}
Selected metric
tonalVar <- "dTonalECMAAvgMaxLR"
# Fluctuation strength
iVars <- c(diffVar, "dFluctSHM10ExBin", "dFluctSHM05ExBin", "dFluctOV10ExMaxLR", "dFluctOV05ExMaxLR")
dVar <- "NoticedPcFilt"
seeds <- c(418657, 84, 1630, 18659, 3687)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 1501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutFluct <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutFluct$OOB_RMSE
[1] 11.38464
resultsOutFluct$OOB_MAE
[1] 8.885649
resultsOutFluct$Rsquared
[1] 0.8846925
Train multiple seeds model
resultsOutFluct <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutFluct$OOB_RMSE
[1] 11.39714
resultsOutFluct$OOB_MAE
[1] 8.868102
resultsOutFluct$Rsquared
[1] 0.8843449
# store results
resNoticeFit['Diff fluct', 'RMSE'] <- resultsOutFluct$OOB_RMSE
resNoticeFit['Diff fluct', 'MAE'] <- resultsOutFluct$OOB_MAE
resNoticeFit['Diff fluct', 'Rsquared'] <- resultsOutFluct$Rsquared
resNoticePermImp$DiffFluct <- resultsOutFluct$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutFluct.conimp <- arrange(resultsOutFluct$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutFluct.conimp) + geom_col(aes(x=factor(rownames(resultsOutFluct.conimp), levels=rownames(resultsOutFluct.conimp)), y=CondPermImp), fill=mycolours[4], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("dFluctuation strength") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtANoticedFluctConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtANoticeFluctConPermimp.svg")
ggsave(filename="PtANoticedFluctConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeFluctConPermimp.pdf")
}
Selected metric
fluctVar <- "dFluctSHM05ExBin"
# Roughness
iVars <- c(diffVar, "dRoughECMA10ExBin", "dRoughECMA05ExBin", "dRoughFZ10ExMaxLR", "dRoughFZ05ExMaxLR")
dVar <- "NoticedPcFilt"
seeds <- c(69851, 85109, 410986, 1563, 896)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 2501
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutRough <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutRough$OOB_RMSE
[1] 12.37658
resultsOutRough$OOB_MAE
[1] 9.87597
resultsOutRough$Rsquared
[1] 0.8723291
Train multiple seeds model
resultsOutRough <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutRough$OOB_RMSE
[1] 12.37394
resultsOutRough$OOB_MAE
[1] 9.849881
resultsOutRough$Rsquared
[1] 0.8718619
# store results
resNoticeFit['Diff rough', 'RMSE'] <- resultsOutRough$OOB_RMSE
resNoticeFit['Diff rough', 'MAE'] <- resultsOutRough$OOB_MAE
resNoticeFit['Diff rough', 'Rsquared'] <- resultsOutRough$Rsquared
resNoticePermImp$DiffRough <- resultsOutRough$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutRough.conimp <- arrange(resultsOutRough$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutRough.conimp) + geom_col(aes(x=factor(rownames(resultsOutRough.conimp), levels=rownames(resultsOutRough.conimp)), y=CondPermImp), fill=mycolours[5], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("dRoughness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtANoticedRoughConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtANoticedRoughConPermimp.svg")
ggsave(filename="PtANoticedRoughConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtANoticedRoughConPermimp.pdf")
}
Selected metric
roughVar <- "dRoughECMA05ExBin"
# Impulsiveness
iVars <- c(diffVar, "dImpulsSHMAvgMaxLR", "dImpulsSHM05ExMaxLR", "dImpulsSHMPowAvgMaxLR")
dVar <- "NoticedPcFilt"
seeds <- c(418659, 7805, 38475, 65834, 1653)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 251
mtry <- as.integer(length(iVars)/1.25)
Train preliminary model
nperm <- 5
resultsOutImpuls <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutImpuls$OOB_RMSE
[1] 11.85234
resultsOutImpuls$OOB_MAE
[1] 9.182759
resultsOutImpuls$Rsquared
[1] 0.8773988
Train multiple seeds model
resultsOutImpuls <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutImpuls$OOB_RMSE
[1] 11.92896
resultsOutImpuls$OOB_MAE
[1] 9.255452
resultsOutImpuls$Rsquared
[1] 0.8757174
# store results
resNoticeFit['Diff impuls', 'RMSE'] <- resultsOutImpuls$OOB_RMSE
resNoticeFit['Diff impuls', 'MAE'] <- resultsOutImpuls$OOB_MAE
resNoticeFit['Diff impuls', 'Rsquared'] <- resultsOutImpuls$Rsquared
resNoticePermImp$DiffImpuls <- resultsOutImpuls$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutImpuls.conimp <- arrange(resultsOutImpuls$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutImpuls.conimp) + geom_col(aes(x=factor(rownames(resultsOutImpuls.conimp), levels=rownames(resultsOutImpuls.conimp)), y=CondPermImp), fill=mycolours[6], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + ggtitle("dImpulsiveness") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip()
pBar
if (saveplots){
ggsave(filename="PtANoticedImpulsConPermimp.svg", width=8, height=2, path=file.path(outFigPath, "svg"))
unlink("PtANoticedImpulsConPermimp.svg")
ggsave(filename="PtANoticedImpulsConPermimp.pdf", width=8, height=2, path=file.path(outFigPath, "pdf"))
unlink("PtANoticedImpulsConPermimp.pdf")
}
Selected metric
impulsVar <- "dImpulsSHM05ExMaxLR"
Now the highest importance dSQMs are ranked against each other, controlling for loudness difference.
iVars <- c(diffVar, dSharpVar, dTonLdVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "NoticedPcFilt"
seeds <- c(98465, 54163, 6541, 36485, 849675)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 501
mtry <- as.integer(length(iVars)/1.75)
Train preliminary model
nperm <- 5
resultsOutSQMs1 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs1$OOB_RMSE
[1] 10.79004
resultsOutSQMs1$OOB_MAE
[1] 8.514653
resultsOutSQMs1$Rsquared
[1] 0.9054054
Train multiple seeds model
resultsOutSQMs1 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs1$OOB_RMSE
[1] 10.5503
resultsOutSQMs1$OOB_MAE
[1] 8.348676
resultsOutSQMs1$Rsquared
[1] 0.9084059
# store results
resNoticeFit['Diff SQMs inc tonal loud', 'RMSE'] <- resultsOutSQMs1$OOB_RMSE
resNoticeFit['Diff SQMs inc tonal loud', 'MAE'] <- resultsOutSQMs1$OOB_MAE
resNoticeFit['Diff SQMs inc tonal loud', 'Rsquared'] <- resultsOutSQMs1$Rsquared
resNoticePermImp$DiffSQMs1 <- resultsOutSQMs1$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSQMs1.conimp <- arrange(resultsOutSQMs1$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSQMs1.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs1.conimp), levels=rownames(resultsOutSQMs1.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 600))
pBar
if (saveplots){
ggsave(filename="PtANoticeDiffSQMsTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
unlink("PtANoticeDiffSQMsTonLdConPermimp.svg")
ggsave(filename="PtANoticeDiffSQMsTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeDiffSQMsTonLdConPermimp.pdf")
}
iVars <- c(diffVar, dSharpVar, dTonalVar, dFluctVar, dRoughVar, dImpulsVar)
dVar <- "NoticedPcFilt"
seeds <- c(49865, 7852, 845961, 410583, 36748)
p <- mtryTune(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seed=seeds[1],
ntrees=ntrees, minsplit=minsplit, minbucket=minbucket)
p
Selected hyperparameters
ntree <- 4001
mtry <- as.integer(length(iVars)/1.75)
Train preliminary model
nperm <- 5
resultsOutSQMs2 <- crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds[1:2], ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs2$OOB_RMSE
[1] 11.8702
resultsOutSQMs2$OOB_MAE
[1] 9.243409
resultsOutSQMs2$Rsquared
[1] 0.8859168
Train multiple seeds model
resultsOutSQMs2 <- multi_crfReg(dataIn=stimNoticeDataANum, iVars=iVars, dVar=dVar, seeds=seeds, ntree=ntree, mtry=mtry, permImpCondThres=permImpCondThres, nperm=nperm, minsplit=minsplit, minbucket=minbucket)
# print model prediction results
resultsOutSQMs2$OOB_RMSE
[1] 11.81577
resultsOutSQMs2$OOB_MAE
[1] 9.205458
resultsOutSQMs2$Rsquared
[1] 0.8863226
# store results
resNoticeFit['Diff SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_RMSE
resNoticeFit['Diff SQMs no tonal loud', 'RMSE'] <- resultsOutSQMs2$OOB_MAE
resNoticeFit['Diff SQMs no tonal loud', 'Rsquared'] <- resultsOutSQMs2$Rsquared
resNoticePermImp$DiffSQMs2 <- resultsOutSQMs2$conditional_permimp
par(mai=c(0,3,0,0))
# plot conditional importance
resultsOutSQMs2.conimp <- arrange(resultsOutSQMs2$conditional_permimp, desc(row_number()))
pBar <- ggplot(resultsOutSQMs2.conimp) + geom_col(aes(x=factor(rownames(resultsOutSQMs2.conimp), levels=rownames(resultsOutSQMs2.conimp)), y=CondPermImp), fill=mycolours[7], width=0.5) + labs(x="Variable", y="Conditional variable permutation importance (% UAS noticed)") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2)) + coord_flip(ylim=c(0, 600))
pBar
if (saveplots){
ggsave(filename="PtANoticeDiffSQMsNoTonLdConPermimp.svg", width=8, height=2.4, path=file.path(outFigPath, "svg"))
unlink("PtANoticeDiffSQMsNoTonLdConPermimp.svg")
ggsave(filename="PtANoticeDiffSQMsNoTonLdConPermimp.pdf", width=8, height=2.4, path=file.path(outFigPath, "pdf"))
unlink("PtANoticeDiffSQMsNoTonLdConPermimp.pdf")
}
if (savedata){
utils::write.csv(resNoticeFit, paste(outDataPath, "\\ptACRFNoticeOOBFit.csv", sep=""))
ii <- 0
temp = list()
for (res in resNoticePermImp){
ii <- ii + 1
temp[[ii]] <- as.data.frame(resNoticePermImp[ii])
names(temp[[ii]]) <- names(resNoticePermImp[ii])
}
openxlsx::write.xlsx(temp, paste(outDataPath, "\\ptACRFNoticeConPermimp.xlsx",
sep=""),
rowNames=TRUE)
}
Summary of results for Part A
# combine the annoyance perm importance results
# convert each result to a tibble with rownames added to a column, renaming the data column to 'dAnnoy' etc.
resdAnnoyMnAbsPermImpTbl <- as.data.frame(resdAnnoyMnPermImp$AbsSQMs1/max(resdAnnoyMnPermImp$AbsSQMs1)) |>
tibble::rownames_to_column(var='Variable')
colnames(resdAnnoyMnAbsPermImpTbl)[2] <- "dAnnoy"
resdHiAnnoyAbsPermImpTbl <- as.data.frame(resdHiAnnoyPermImp$AbsSQMs1/max(resdHiAnnoyPermImp$AbsSQMs1)) |>
tibble::rownames_to_column(var='Variable')
colnames(resdHiAnnoyAbsPermImpTbl)[2] <- "dHiAnnoy"
resNoticeAbsPermImpTbl <- as.data.frame(resNoticePermImp$AbsSQMs1/max(resNoticePermImp$AbsSQMs1)) |>
tibble::rownames_to_column(var='Variable')
colnames(resNoticeAbsPermImpTbl)[2] <- "dNotice"
# merge the dataframes
resAbsPermImpTbl <- list(resdAnnoyMnAbsPermImpTbl, resdHiAnnoyAbsPermImpTbl, resNoticeAbsPermImpTbl) |>
purrr::reduce(merge, by = c('Variable'), all = T)
# rename the columns
colnames(resAbsPermImpTbl)[2:4] <- c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed")
resAbsPermImpTbl[is.na(resAbsPermImpTbl)] <- 0
resAbs <- tidyr::pivot_longer(resAbsPermImpTbl, cols=-Variable, names_to="Outcome", values_to="Imp")
# reorder res tibble, descending by the variable Imp grouped sum and create column with new group order as a factor
resAbs <- resAbs |> mutate(Variable_sum = sum(Imp), .by=Variable) |> arrange(desc(Variable_sum)) |> group_by(Variable_sum, Variable) |>
mutate(Order = cur_group_id()) |> mutate(Order = as.factor(Order)) |> arrange(desc(Order))
# Reorder outcome levels
resAbs$Outcome <- factor(resAbs$Outcome, levels=c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed"))
# plot res as horizontal bar chart, with Imp as y axis, Variable as x axis, Outcome as fill, and Variable_sum as order, relabel x axis with Variable names
pBar <- ggplot(resAbs) + geom_col(aes(fill=Outcome, y=Imp, x=Order), colour='grey35', width=0.75, show.legend=TRUE) + labs(x="Variable", y="Normalised variable permutation importance") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2), legend.position = "right") + coord_flip(ylim=c(-0.1, 2)) + scale_fill_manual(values=mycolours) + scale_x_discrete(labels=unique(rev(resAbs$Variable)))
pBar
if (saveplots){
ggsave(filename="PtAcrfAbsSQMsSummary.svg", width=8, height=4, path=file.path(outFigPath, "svg"))
unlink("PtAcrfAbsSQMsSummary.svg")
ggsave(filename="PtAcrfAbsSQMsSummary.pdf", width=8, height=4, path=file.path(outFigPath, "pdf"))
unlink("PtAcrfAbsSQMsSummary.pdf")
}
# combine the annoyance perm importance results
# convert each result to a tibble with rownames added to a column, renaming the data column to 'dAnnoy' etc.
resdAnnoyMnDiffPermImpTbl <- as.data.frame(resdAnnoyMnPermImp$DiffSQMs1/max(resdAnnoyMnPermImp$DiffSQMs1)) |>
tibble::rownames_to_column(var='Variable')
colnames(resdAnnoyMnDiffPermImpTbl)[2] <- "dAnnoy"
resdHiAnnoyDiffPermImpTbl <- as.data.frame(resdHiAnnoyPermImp$DiffSQMs1/max(resdHiAnnoyPermImp$DiffSQMs1)) |>
tibble::rownames_to_column(var='Variable')
colnames(resdHiAnnoyDiffPermImpTbl)[2] <- "dHiAnnoy"
resNoticeDiffPermImpTbl <- as.data.frame(resNoticePermImp$DiffSQMs1/max(resNoticePermImp$DiffSQMs1)) |>
tibble::rownames_to_column(var='Variable')
colnames(resNoticeDiffPermImpTbl)[2] <- "dNotice"
# merge the dataframes
resDiffPermImpTbl <- list(resdAnnoyMnDiffPermImpTbl, resdHiAnnoyDiffPermImpTbl, resNoticeDiffPermImpTbl) |>
purrr::reduce(merge, by = c('Variable'), all = T)
# rename the columns
colnames(resDiffPermImpTbl)[2:4] <- c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed")
resDiffPermImpTbl[is.na(resDiffPermImpTbl)] <- 0
resDiff <- tidyr::pivot_longer(resDiffPermImpTbl, cols=-Variable, names_to="Outcome", values_to="Imp")
# reorder res tibble, descending by the variable Imp grouped sum and create column with new group order as a factor
resDiff <- resDiff |> mutate(Variable_sum = sum(Imp), .by=Variable) |> arrange(desc(Variable_sum)) |> group_by(Variable_sum, Variable) |>
mutate(Order = cur_group_id()) |> mutate(Order = as.factor(Order)) |> arrange(desc(Order))
# Reorder outcome levels
resDiff$Outcome <- factor(resDiff$Outcome, levels=c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed"))
# plot res as horizontal bar chart, with Imp as y axis, Variable as x axis, Outcome as fill, and Variable_sum as order, relabel x axis with Variable names
pBar <- ggplot(resDiff) + geom_col(aes(fill=Outcome, y=Imp, x=Order), colour='grey35', width=0.75, show.legend=TRUE) + labs(x="Variable", y="Normalised variable permutation importance") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2), legend.position = "right") + coord_flip(ylim=c(-0.1, 2)) + scale_fill_manual(values=mycolours) + scale_x_discrete(labels=unique(rev(resDiff$Variable)))
pBar
if (saveplots){
ggsave(filename="PtAcrfDiffSQMsSummary.svg", width=8, height=4, path=file.path(outFigPath, "svg"))
unlink("PtAcrfDiffSQMsSummary.svg")
ggsave(filename="PtAcrfDiffSQMsSummary.pdf", width=8, height=4, path=file.path(outFigPath, "pdf"))
unlink("PtAcrfDiffSQMsSummary.pdf")
}
# combine the annoyance perm importance results
# convert each result to a tibble with rownames added to a column, renaming the data column to 'dAnnoy' etc.
resdAnnoyMnAbsPermImpNoTonLdTbl <- as.data.frame(resdAnnoyMnPermImp$AbsSQMs2/max(resdAnnoyMnPermImp$AbsSQMs2)) |>
tibble::rownames_to_column(var='Variable')
colnames(resdAnnoyMnAbsPermImpNoTonLdTbl)[2] <- "dAnnoy"
resdHiAnnoyAbsPermImpNoTonLdTbl <- as.data.frame(resdHiAnnoyPermImp$AbsSQMs2/max(resdHiAnnoyPermImp$AbsSQMs2)) |>
tibble::rownames_to_column(var='Variable')
colnames(resdHiAnnoyAbsPermImpNoTonLdTbl)[2] <- "dHiAnnoy"
resNoticeAbsPermImpNoTonLdTbl <- as.data.frame(resNoticePermImp$AbsSQMs2/max(resNoticePermImp$AbsSQMs2)) |>
tibble::rownames_to_column(var='Variable')
colnames(resNoticeAbsPermImpNoTonLdTbl)[2] <- "dNotice"
# merge the dataframes
resAbsNoTonLdPermImpNoTonLdTbl <- list(resdAnnoyMnAbsPermImpNoTonLdTbl, resdHiAnnoyAbsPermImpNoTonLdTbl, resNoticeAbsPermImpNoTonLdTbl) |>
purrr::reduce(merge, by = c('Variable'), all = T)
# rename the columns
colnames(resAbsNoTonLdPermImpNoTonLdTbl)[2:4] <- c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed")
resAbsNoTonLdPermImpNoTonLdTbl[is.na(resAbsNoTonLdPermImpNoTonLdTbl)] <- 0
resAbsNoTonLd <- tidyr::pivot_longer(resAbsNoTonLdPermImpNoTonLdTbl, cols=-Variable, names_to="Outcome", values_to="Imp")
# reorder res tibble, descending by the variable Imp grouped sum and create column with new group order as a factor
resAbsNoTonLd <- resAbsNoTonLd |> mutate(Variable_sum = sum(Imp), .by=Variable) |> arrange(desc(Variable_sum)) |> group_by(Variable_sum, Variable) |>
mutate(Order = cur_group_id()) |> mutate(Order = as.factor(Order)) |> arrange(desc(Order))
# Reorder outcome levels
resAbsNoTonLd$Outcome <- factor(resAbsNoTonLd$Outcome, levels=c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed"))
# plot res as horizontal bar chart, with Imp as y axis, Variable as x axis, Outcome as fill, and Variable_sum as order, relabel x axis with Variable names
pBar <- ggplot(resAbsNoTonLd) + geom_col(aes(fill=Outcome, y=Imp, x=Order), colour='grey35', width=0.75, show.legend=TRUE) + labs(x="Variable", y="Normalised variable permutation importance") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2), legend.position = "right") + coord_flip(ylim=c(-0.1, 2)) + scale_fill_manual(values=mycolours) + scale_x_discrete(labels=unique(rev(resAbsNoTonLd$Variable)))
pBar
if (saveplots){
ggsave(filename="PtAcrfAbsSQMsNoTonLdSummary.svg", width=8, height=3.7, path=file.path(outFigPath, "svg"))
unlink("PtAcrfAbsSQMsNoTonLdSummary.svg")
ggsave(filename="PtAcrfAbsSQMsNoTonLdSummary.pdf", width=8, height=3.7, path=file.path(outFigPath, "pdf"))
unlink("PtAcrfAbsSQMsNoTonLdSummary.pdf")
}
# combine the annoyance perm importance results
# convert each result to a tibble with rownames added to a column, renaming the data column to 'dAnnoy' etc.
resdAnnoyMnDiffPermImpNoTonLdTbl <- as.data.frame(resdAnnoyMnPermImp$DiffSQMs2/max(resdAnnoyMnPermImp$DiffSQMs2)) |>
tibble::rownames_to_column(var='Variable')
colnames(resdAnnoyMnDiffPermImpNoTonLdTbl)[2] <- "dAnnoy"
resdHiAnnoyDiffPermImpNoTonLdTbl <- as.data.frame(resdHiAnnoyPermImp$DiffSQMs2/max(resdHiAnnoyPermImp$DiffSQMs2)) |>
tibble::rownames_to_column(var='Variable')
colnames(resdHiAnnoyDiffPermImpNoTonLdTbl)[2] <- "dHiAnnoy"
resNoticeDiffPermImpNoTonLdTbl <- as.data.frame(resNoticePermImp$DiffSQMs2/max(resNoticePermImp$DiffSQMs2)) |>
tibble::rownames_to_column(var='Variable')
colnames(resNoticeDiffPermImpNoTonLdTbl)[2] <- "dNotice"
# merge the dataframes
resDiffNoTonLdPermImpNoTonLdTbl <- list(resdAnnoyMnDiffPermImpNoTonLdTbl, resdHiAnnoyDiffPermImpNoTonLdTbl, resNoticeDiffPermImpNoTonLdTbl) |>
purrr::reduce(merge, by = c('Variable'), all = T)
# rename the columns
colnames(resDiffNoTonLdPermImpNoTonLdTbl)[2:4] <- c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed")
resDiffNoTonLdPermImpNoTonLdTbl[is.na(resDiffNoTonLdPermImpNoTonLdTbl)] <- 0
resDiffNoTonLd <- tidyr::pivot_longer(resDiffNoTonLdPermImpNoTonLdTbl, cols=-Variable, names_to="Outcome", values_to="Imp")
# reorder res tibble, descending by the variable Imp grouped sum and create column with new group order as a factor
resDiffNoTonLd <- resDiffNoTonLd |> mutate(Variable_sum = sum(Imp), .by=Variable) |> arrange(desc(Variable_sum)) |> group_by(Variable_sum, Variable) |>
mutate(Order = cur_group_id()) |> mutate(Order = as.factor(Order)) |> arrange(desc(Order))
# Reorder outcome levels
resDiffNoTonLd$Outcome <- factor(resDiffNoTonLd$Outcome, levels=c("Mean change in annoyance", "%HA [p(HA | 0)]", "% UAS noticed"))
# plot res as horizontal bar chart, with Imp as y axis, Variable as x axis, Outcome as fill, and Variable_sum as order, relabel x axis with Variable names
pBar <- ggplot(resDiffNoTonLd) + geom_col(aes(fill=Outcome, y=Imp, x=Order), colour='grey35', width=0.75, show.legend=TRUE) + labs(x="Variable", y="Normalised variable permutation importance") + theme(text = element_text(family = "serif"), panel.grid=element_line(color = rgb(235, 235, 235, 100, maxColorValue = 255), linewidth = 0.25, linetype = 2), legend.position = "right") + coord_flip(ylim=c(-0.1, 2)) + scale_fill_manual(values=mycolours) + scale_x_discrete(labels=unique(rev(resDiffNoTonLd$Variable)))
pBar
if (saveplots){
ggsave(filename="PtAcrfDiffSQMsNoTonLdSummary.svg", width=8, height=4, path=file.path(outFigPath, "svg"))
unlink("PtAcrfDiffSQMsNoTonLdSummary.svg")
ggsave(filename="PtAcrfDiffSQMsNoTonLdSummary.pdf", width=8, height=4, path=file.path(outFigPath, "pdf"))
unlink("PtAcrfDiffSQMsNoTonLdSummary.pdf")
}
# Make a list of the summary results
resSummary <- list(resAbs, resDiff, resAbsNoTonLd, resDiffNoTonLd)
# Save the results
if (savedata){
ii <- 0
temp = list()
for (res in resSummary){
ii <- ii + 1
temp[[ii]] <- data.frame(resSummary[ii])
}
openxlsx::write.xlsx(temp, paste(outDataPath, "\\ptACRFSummary.xlsx",
sep=""),
rowNames=TRUE)
}